{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solver tolerances\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import pybamm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The goal of each solver is the solve the problem as fast as possible while maintaining a certain level of accuracy. The accuracy is controlled by the tolerances, one for the relative residual and one for the absolute residual.\n",
    "\n",
    "At each iteration, the solver computes an estimate of the current error $\\mathbf{e_k} = \\mathbf{y_k} - \\mathbf{y}$, where $\\mathbf{y_k}$ is the current approximation and $\\mathbf{y}$ is the exact solution. It accepts the current approximation if the error is small enough according to the set relative tolerance $\\epsilon_r$ and absolute tolerance $\\epsilon_a$, using the following criteria:\n",
    "\n",
    "$$\n",
    "\\sqrt{\\frac{1}{n} \\sum_{i=1}^n \\left( \\frac{e_{k,i}}{\\epsilon_r |y_i| + \\epsilon_a} \\right)^2 } \\leq 1\n",
    "$$\n",
    "\n",
    "The relative tolerance is used to control the error relative to the current magnitude of the solution, while the absolute tolerance is used to control the error in absolute terms (i.e. when the solution is close to zero). Note that these errors are computed in relation to the solution for the state variables of the equations, and not the output variables that might be computed by PyBaMM after the solution is found.\n",
    "\n",
    "## Default tolerances\n",
    "\n",
    "The default tolerances for the solvers are always given in the documentation of each solver. For example, the default tolerances for the IDAKLU solver can be seen [here](https://docs.pybamm.org/en/stable/source/api/solvers/idaklu_solver.html) as the default values for the `rtol` and `atol` arguments to the solver `__init__` method.\n",
    "\n",
    "## Setting tolerances\n",
    "\n",
    "Generally, a smaller value for the tolerances will lead to a more accurate solution, but will also require more computational effort. For examples, here are the timings for the default SPM and DFN models with different tolerances ranging from $10^{-2}$ to $10^{-6}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHtCAYAAAAEKx0CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABIqUlEQVR4nO3dd3hUZfrG8XuSkEYKNUCogdBrEFQEpGpAREXBVUEQ1BUFBWnKT1GsCBpxV1yK7oKFIkhZpahZMCAaRZBQpEOE0BIUSIE0Juf3R2R0SCGEvMwQvp/rmusy77zvmedwJA/3zDlnbJZlWQIAAAAAACXOw9UFAAAAAABQWhG6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDaDItm3bpr59+6p27dry9fVV9erVdcstt+jdd991zKlTp45sNpvjERISoo4dO2rp0qVO2+rcubNsNpvq16+f72tFR0c7tvHZZ58Z3S8AAEqTOXPmOPViX19fhYaGKjIyUv/85z+VmprqNH/ixIlO8//6mDFjhmPe+bGoqKgCX3Pjxo3G9w+42ni5ugAAV4fvv/9eXbp0Ua1atfToo4+qatWqSkhI0A8//KB//OMfevLJJx1zW7VqpdGjR0uSjh49qpkzZ+ruu+/W9OnTNXToUMc8X19f7du3Txs2bND111/v9Hpz586Vr6+vMjIyrswOAgBQyrz88ssKCwtTdna2jh8/rpiYGI0cOVJvv/22Pv/8c7Vo0cJp/vTp0xUQEOA0dsMNN+TZ7ptvvqnHH39c/v7+RusHSgtCN4Aiee211xQcHKyffvpJ5cqVc3ouKSnJ6efq1atrwIABjp8HDhyo8PBwTZ061Sl016tXT+fOndP8+fOdQndGRoaWLl2qXr16afHixWZ2CACAUq5nz55q06aN4+fx48drzZo1uv3223XHHXdo586d8vPzczzft29fVapUqdBttmrVSnFxcZoxY4ZGjRplrHagNOH0cgBFsn//fjVt2jRP4JakkJCQQtdWrVpVjRs3Vnx8fJ7n7r//fn366afKyclxjH3xxRc6e/as7r333suuGwAA/Klr166aMGGCDh48qE8++eSS17dv315du3bVlClTlJ6ebqBCoPQhdAMoktq1a2vTpk3avn37Ja/Nzs5WQkKCKlasmOe5Bx54QMeOHVNMTIxjbN68eerWrdtFwzwAALh0Dz74oCTp66+/dho/efKkfvvtN8fj1KlT+a6fOHGiEhMTNX36dOO1AqUBoRtAkYwZM0Znz55Vq1atdNNNN+mZZ57R119/rezs7Dxzs7OzHQ1769atGjhwoBITE9WvX788c+vXr682bdpo3rx5kqTTp09r5cqVeuCBB4zvEwAA16IaNWooODhY+/fvdxpv2LChKleu7HhERETku75jx47q0qWL3nzzTT7tBoqA0A2gSG655RbFxsbqjjvu0JYtWzRlyhRFRkaqevXq+vzzz53mfv31146G3bJlSy1atEgPPvigJk+enO+2H3jgAS1ZskRZWVn67LPP5OnpqT59+lyJ3QIA4JoUEBCQ5y7mixcvVnR0tOMxd+7cAtdPnDhRx48fd7q7OYD8cSM1AEXWtm1bRzjesmWLli5dqqlTp6pv376Ki4tTkyZNJOXe6fTVV1+VzWaTv7+/GjdunO+14Ofdd999GjNmjFatWqW5c+fq9ttvV2Bg4BXaKwAArj1paWl5LuO6+eabL3ojtb/O7dKli6ZMmeJ0k1QAeRG6AVwyb29vtW3bVm3btlWDBg00ePBgLVq0SC+++KIkqVKlSurevXuRt1etWjV17txZUVFR+u6777hjOQAABh0+fFjJyckKDw+/rO28+OKL6ty5s2bOnFnom+vAtY7TywFclvNfRXLs2LHL2s4DDzygb7/9VkFBQbrttttKojQAAJCPjz/+WJIUGRl5Wdvp1KmTOnfurMmTJ3NtN1AIPukGUCTffPONOnfuLJvN5jS+cuVKSbk3X7kcffv2VUJCgho2bChvb+/L2hYAAMjfmjVr9MorrygsLEz9+/e/7O1NnDhRnTt31qxZs0qgOqB0InQDKJInn3xSZ8+eVZ8+fdSoUSNlZWXp+++/16effqo6depo8ODBl7X94OBgTZw4sWSKBQAAWrVqlXbt2qVz584pMTFRa9asUXR0tGrXrq3PP/9cvr6+l/0anTp1UqdOnbR27doSqBgonQjdAIrkrbfe0qJFi7Ry5UrNmjVLWVlZqlWrlp544gk9//zzXMsFAICbeeGFFyTl3oulQoUKat68ud555x0NHjy4RG9YOnHiRHXp0qXEtgeUNjbLsixXFwEAAAAAQGnEjdQAAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABD+MqwYsrJydHRo0cVGBgom83m6nIAAKWIZVlKTU1VaGioPDx4f/xy0bMBACYUtV8Tuovp6NGjqlmzpqvLAACUYgkJCapRo4ary7jq0bMBACZdrF8TuospMDBQktRBt8lLZVxcDZx4eLq6AuTD5slxcVdWdparS8AFzilb67XS0Wtwec7/Od7s00deNnq2W+FMDreUczbd1SUAV4Wi9mtCdzGdPz3NS2Vo4O7GRrhzRzaOi9uybJarS8CF/jgknApdMhw921ZGXjZvF1cDJzZCtzvKsZ1zdQnA1aGI/ZrfdAAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAY4uXqAq51Zby9NOjlv6n7gJsVUD5AB7Ye1JwJC/Tz/7YWaX3F0Ap6/O1Buu7WlrJ52LTlm180fdQcHY9PyjO3x5Cu6je6t6qGhehEwu9a+u5K/Xfal8XeZnTOonxr+mD8XH06eVmR6gcA4GpRxttLAyfco24PtFdAubKK356gD1/6TD+v2V6k9RVDy2vo5P5q3a2ZbB4e2rpuh2aMm6vjv55wmnf7o93UslMTNWpbVyE1K+nrj79V1GOz8myvQtVg3fVEpBq1raf6rcPkH+insZGvaeu3u5zmValVSR/tmlpgXatmf6N3hv2nSPsAALh0hG4XGzt7mDr2vVFL/rFSR/Ye062DOuu1FeM1putL+uW7XYWu9S3rq7fWvKiywf6aP2mJzmXbdc/I2xUV85KGRoxV6sk0x9xef++ukTMe07rPftBnU5ereYfGGv7Ph+Xr76NPp/y3WNuUpE1fb1H0x2udxvZtji+BPxkAANzL6Fl/V8c+bbV02lc6sj9Rtw7oqFeWjta4HpP0S+yeQtf6lvXRlFXjVTbIXwve/ELnsu26+8lIvfX1c3r8xued+mu/Ub3kH+Cr3ZsOqELVcgVus0b9avrbmN46vPeYfv3lsJrcWD/fead/S9XkIdPzjLe5pYW63d9em/5XtDcNAADFQ+h2oYZtw9Xl/g6aOfYjfRb1hSQp+qO1en/b23p08gCN7PB8oevveOJW1WgQqmHXP6s9G/dLkn5atVnvb3tb/Ub31n+emy9J8vb11uBX79cPyzfplXujJEmrPlgtm4dN/Z/vqxWz/qe002cuaZvnHd57TKvnfltyfygAALihhm3qqsu97fT++Pn67B8rJUn/m7teszZO0iOv3aenu75c6Pref++uGvWr6cmOL2jPptw3pzd+vUUzN05S3xE9NfvFP88eG3vra0pK+F2StCzp/QK3uXfzr+pbfahST51Rh7vaFhi6M89mas2C7/OM3zqgo84kn9UPKzcXvvMAgMvCNd0u1LHvjbKfs2vlrP85xrIzs/Xlf1ar6U0NVblGxcLX39NOuzbsc4RjSUrYfVSbV2/Tzf1ucoy16tJUwZWC9MX0r5zWf/6vr+QX4KsberW+5G3+lbevt8r4lCnaTv/BL8BXj099SB8feE8r0udp4fEP9MZXExQeEXZJ2wEA4ErocNf1uT37P2scY9mZ2fryw7VqcmN9Va5eodD1Hfu01e6N+x2BW5IS9hzT5m9+0c133+A093zgvpj0tAylnjpzCXvxpwpVg9WiUxN99/lGZWdmX3T+HUNv0ayNk/Tf3z7QZ0dm6N31L6nLve2K9doAcK1xaehet26devfurdDQUNlsNi1btqxEthsTE6PWrVvLx8dH4eHhmjNnjtPzkyZNUtu2bRUYGKiQkBDddddd2r17d4m89qUIbxWmw3uO6WxqutP47g37JEn1WtUpcK3NZlPdFrW0Z9P+PM/t+mmfqodXlV+Ab+7r/BFk/xqkJWnvpgOy23Mcz1/KNs+7dVBnfZH2sVamz9MH26eqy/0dLrLXuUZM/7tuH3qr1i/5Ue8O+0CfRX2urPQs1WpcvUjrAQBX1jXfs1vW1uG9x3U2NcNpfPcfvbVuy9oFrrXZbAprVlN7fs57+dXuTQcUWq9Knv5qWqe+7eTp6ZHvJ+AX6jm4s4a9PVAHdx3VjHFz9fGrS7R/6yE1bFvPfKEAUAq49PTyM2fOqGXLlhoyZIjuvvvuEtlmfHy8evXqpaFDh2ru3LlavXq1HnnkEVWrVk2RkZGSpLVr12rYsGFq27atzp07p//7v//Trbfeqh07dqhs2bIlUkdRVKhWTiePncoz/vux05Jyb7hSkMAKAfL29c53/fmxiqEVdHjPUVWoVl72c3adPpHiNO9c9jml/J6qitXKX/I2JemX73Zp7aJYHY9PUsXQ8rrjiR76v7kjVDbYX8tnfF3ovt/Qq7VWfbBaM8d85Bhb+Obnha4BALjONd+zq5bTyeOn84yfH6tYrVyBawMrlM3trxdZf3jv8csvtIi63tdOvx87pbiYHRede32PVvr1l8N6bcC7V6AyACh9XBq6e/bsqZ49exb4fGZmpp577jnNnz9fp0+fVrNmzTR58mR17ty5wDUzZsxQWFiYoqJyr11u3Lix1q9fr6lTpzoa+JdfOt+xe86cOQoJCdGmTZt08803F1hLZmam4+eUlJR8510KHz/vfE/pysrIkiR5+3kXulZSAeuzneZ4+3krO+tcvtvJzsh2vM6lbFOSRnac4DTny/98o39tnKwhr92vr+fEOPYjP2mnz6jR9eGqWK28fs8n5AMA3Mu13rO9/crk20sd/dG3kJ7te76/Fry+sJ5f0qqHV1WD1nW1+J+rZFnWReennT6rStXLq8F1YU6nxwMAisatr+kePny4YmNjtWDBAm3dulX9+vVTjx49tHfv3gLXxMbGqnv37k5jkZGRio2NLXBNcnKyJKlChYKvx5o0aZKCg4Mdj5o1a17i3uSVmZ6V77XQ3n8056z0gkNr5h/P5b++jNOcrPQslfHO//2VMr5lHK9zKdvMz7nsc/rve6sUWD5A9a+rW+A8SXr/mU9Up1ktzT00Q+/+MEkPvthPVcNCCl0DAHBfpb1nZ6Vn59tLHf2xkDeazz9Xxqfg9YX1/JLW9b7ce7Ss+fTip5ZL0sK3lyv9TKbe/fZl/Wfrmxo2dVCBN20DAOTltqH70KFDmj17thYtWqSOHTuqXr16GjNmjDp06KDZs2cXuO748eOqUqWK01iVKlWUkpKi9PT0PPNzcnI0cuRItW/fXs2aNStwu+PHj1dycrLjkZCQUPyd+8PJY6dVoVreU8jPn6L2+9GCPwFOPZmmrIysfNefH/v96Mk/XueUPL08Va5ykNM8rzJeCqoY6Pik+VK2WZATf9z8JahCQKHz1i2K1cB6w/TeU//W70dPqt+YO/TB9qlq26NVoesAAO7nmujZx0/n+/Vd58fOXxqWn9STZ3L7azHXl7Qu97ZTwu6j2rf51yLNT9h9VA+3HKfXH5ymX2L3qMNdbTR1zQt68PmSucwAAEo7t/3KsG3btslut6tBgwZO45mZmapYMfeu3gEBfwa7AQMGaMaMGZf8OsOGDdP27du1fv36Quf5+PjIx8fnkrdfmP1bflWrLk3lH+jndDO1Rjfkvnu8P+7XAtdalqX4bYfU4Lq8NzFpfH19Hd1/XOlpuTd72ffHdhq0qacNq/78WpAGberK09PD8TqXss2CVKub+4+nC68fz8/J46f1xfSv9cX0r1WucpD+tWmKHvi/e/TTl3EXXQsAcB/XRM/eelAtOzWWf6Cv083UGv1xM7EDWw4WuNayLMX/clgNWuf9ho5Gberp6IHEi/bXktKwbT1VD6+qD1/+7JLWZZ7N1NrFP2rt4h/lVcZTLywYofvH3aEFb35RpLufA8C1zG1Dd1pamjw9PbVp0yZ5eno6PXe+ccfFxTnGgoJyP8WtWrWqEhMTneYnJiYqKChIfn5+TuPDhw/X8uXLtW7dOtWoUcPAXhTu289ide+YO3Tb37s7vqe7jLeXIh/qop0/7NGJw39+ZUjlmpXk6++thN1H/1y/+Ac98sYANbiurvZsOiBJqtEgVK26NtOiP7YnSXFrtivl91TdPvRWp9Dde2ik0s9k6McVP1/yNoMrBSn5N+dg7Rfgqz4jeun0iRTt/WNtfjw8POQb4KuzKWcdY6dPpOj3o6fyPfUOAODeromevewn9Xu6l24b0tXxPd1lvL1064M3a+eGfTpx5M8zwSrXqJjbs/ccc4ytX7pBD796n+q3DtPeP+5iXqN+VbXq3ESfvbPyiu3H+a/5+ubTgk/hv1BghQClnkxz/Hwu266DO4+qza0t5VXGk9ANABfhtgknIiJCdrtdSUlJ6tixY75zwsPD84y1a9dOK1c6N6/o6Gi1a/fnd0lalqUnn3xSS5cuVUxMjMLCXPPd0Ls27NPahd/r4dcfULmQYB3dd1y3DuykKnUqK+qR6U5zn/lwuFp2bqpbPPo5xj7/11fq+Uh3vbp8vBZFfaFz2efU9+neOpWY7AjxUu6N2ea8sEBPvfeoJnw6Shu/jlOzDo3V/cGb9Z/n5in1VNolb/OOYZFqf+f1il2+UUmHflPFauUVObiLQmpV0uSB7+pcdv43bpMkv0BfzU+YqW8/+0H7t/6q9LQMte7WQo2uD9eM0R+WxB8tAOAKuhZ69u6f9mvd4h81+OV+Cq4cpKMHEnVL/w6qUruS3n78A6e5Yz94TC1vbqxI/wcdY1/MWq2eg7volcWj9dk/VsqebdfdT/XQqaRkLf7nKqf1N9wWobrNa0mSPMt4Kqx5Td3/zJ2SpB9W/Kz47X+eLn9+vPYfX7nZ7YEOanpTQ0nS/Mn/ddquh4dNnfreoB0/7tWx+KQi7/ukL8bpZGKydsTu0amkFNVqGKo7hnbXhi/jrtgn9ABwNXNp6E5LS9O+ffscP8fHxysuLk4VKlRQgwYN1L9/fw0cOFBRUVGKiIjQiRMntHr1arVo0UK9evXKd5tDhw7VtGnTNG7cOA0ZMkRr1qzRwoULtWLFCsecYcOGad68efrvf/+rwMBAHT+e+xUdwcHBed5ZN23yoGl66NB96j7gZgWWL6sDWw9pQu83tO3bnRddm56WoTFdXtTjbz+k/s/dI5uHTVtiftGMUR/m+RT6i+lf61y2XX1H9daNd7TRiYTf9a+nZ2vpP1YWa5u/fLdbTds1VM+HuymoYqAyzmRo94Z9inp4uuK+2V5o3Zlns/T59K/U5paWan/39fLw8NDRfcf1jyfev+hXjQEAXIOeLU15ZKYGvXCPuj3QXoHl/BW/PUEv3PO2tn938e8NT0/L0Nger+mxyQP0wDN3yuZh09Zvd2nmuE+U/Fuq09wOd7bVrQ/++eZF/VZ1VL9VHUnSb0dOOoXuh17s67S2x6BOjv++MHRHdG2mClXKaf7kS/uKzhX//kZd/3aT7n6yp/wCfPTbkVNa9q9ozZ+87JK2AwDXKptVlO+KMCQmJkZdunTJMz5o0CDNmTNH2dnZevXVV/XRRx/pyJEjqlSpkm688Ua99NJLat68eaHbffrpp7Vjxw7VqFFDEyZM0EMPPeR43maz5btu9uzZTvMKk5KSouDgYHXWnfKy5b3bN1zIw/Pic3DF2Tw5Lu7Kyr5yd01G0ZyzshWj/yo5OdlxKrarlYae3dX3XnnZrtxXc6EIPNz2nr7XtJyzZy8+CUCR+7VLQ/fVjNDtxgjdbonQ7b4I3e7HHUP31YzQ7cYI3W6J0A0UTVH7Nb/pAAAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADPFydQFXOw9/P3nYvF1dBv7Cahjm6hKQjxw/ft0ARWWdy5B+/K+ryyh1bAGBsnnQs91JTo0QV5eAfGRV8Xd1CSiAV1q2q0vAX+Scy5BiL96v+aQbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADHFp6F63bp169+6t0NBQ2Ww2LVu2rES2GxMTo9atW8vHx0fh4eGaM2dOgXPfeOMN2Ww2jRw5skReGwCA0oieDQBA8bg0dJ85c0YtW7bUe++9V2LbjI+PV69evdSlSxfFxcVp5MiReuSRR/TVV1/lmfvTTz9p5syZatGiRYm9PgAApRE9GwCA4nFp6O7Zs6deffVV9enTJ9/nMzMzNWbMGFWvXl1ly5bVDTfcoJiYmEK3OWPGDIWFhSkqKkqNGzfW8OHD1bdvX02dOtVpXlpamvr376/3339f5cuXL6ldAgCgVKJnAwBQPG59Tffw4cMVGxurBQsWaOvWrerXr5969OihvXv3FrgmNjZW3bt3dxqLjIxUbGys09iwYcPUq1evPHMLkpmZqZSUFKcHAADIRc8GACB/bhu6Dx06pNmzZ2vRokXq2LGj6tWrpzFjxqhDhw6aPXt2geuOHz+uKlWqOI1VqVJFKSkpSk9PlyQtWLBAP//8syZNmlTkeiZNmqTg4GDHo2bNmsXbMQAAShl6NgAABXPb0L1t2zbZ7XY1aNBAAQEBjsfatWu1f/9+SXIaHzp0aJG2m5CQoBEjRmju3Lny9fUtcj3jx49XcnKy45GQkFCs/QIAoLShZwMAUDAvVxdQkLS0NHl6emrTpk3y9PR0ei4gIECSFBcX5xgLCgqSJFWtWlWJiYlO8xMTExUUFCQ/Pz9t2rRJSUlJat26teN5u92udevWadq0acrMzMzzepLk4+MjHx+fkto9AABKDXo2AAAFc9vQHRERIbvdrqSkJHXs2DHfOeHh4XnG2rVrp5UrVzqNRUdHq127dpKkbt26adu2bU7PDx48WI0aNdIzzzyTb/MGAAAFo2cDAFAwl4butLQ07du3z/FzfHy84uLiVKFCBTVo0ED9+/fXwIEDFRUVpYiICJ04cUKrV69WixYt1KtXr3y3OXToUE2bNk3jxo3TkCFDtGbNGi1cuFArVqyQJAUGBqpZs2ZOa8qWLauKFSvmGQcAALno2QAAFI9Lr+neuHGjIiIiFBERIUkaNWqUIiIi9MILL0iSZs+erYEDB2r06NFq2LCh7rrrLv3000+qVatWgdsMCwvTihUrFB0drZYtWyoqKkoffPCBIiMjr8g+AQBQGtGzAQAoHptlWZari7gapaSkKDg4WF3975OXzdvV5eAvrIZhri4B+cjxc9urWQC3c+5chtb++JqSk5Md1z+j+M737G6VHpaXBz3bneTUCHF1CchHZhV/V5eAAnilZbu6BPzFuXMZWhf76kX7tdvevRwAAAAAgKsdoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDvIo6cevWrUXeaIsWLYpVDAAAuDz0awAA3EuRQ3erVq1ks9lkWVa+z59/zmazyW63l1iBAACg6OjXAAC4lyKH7vj4eJN1AACAEkC/BgDAvRQ5dNeuXdtkHQAAoATQrwEAcC9FDt0X2r9/v9555x3t3LlTktSkSRONGDFC9erVK7HiAADA5aFfAwDgWsW6e/lXX32lJk2aaMOGDWrRooVatGihH3/8UU2bNlV0dHRJ1wgAAIqBfg0AgOsV65PuZ599Vk8//bTeeOONPOPPPPOMbrnllhIpDgAAFB/9GgAA1yvWJ907d+7Uww8/nGd8yJAh2rFjx2UXBQAALh/9GgAA1ytW6K5cubLi4uLyjMfFxSkkJORyawIAACWAfg0AgOsV6/TyRx99VH//+9914MAB3XTTTZKk7777TpMnT9aoUaNKtEB351GhvDw8fFxdBv4itVZZV5eAfCREuroCFMRm2VxdAi6Qky7px8vfDv3amc3LUzYPT1eXgb9IrR/o6hKQj+VRb7u6BBSgXexjri4Bf2E/a5diLz6vWKF7woQJCgwMVFRUlMaPHy9JCg0N1cSJE/XUU08VZ5MAAKCE0a8BAHC9YoVum82mp59+Wk8//bRSU1MlSYGBvFMJAIA7oV8DAOB6xbqmu2vXrjp9+rSk3OZ9voGnpKSoa9euJVYcAAAoPvo1AACuV6zQHRMTo6ysrDzjGRkZ+vbbby+7KAAAcPno1wAAuN4lnV6+detWx3/v2LFDx48fd/xst9v15Zdfqnr16iVXHQAAuGT0awAA3Mclhe5WrVrJZrPJZrPle1qan5+f3n333RIrDgAAXDr6NQAA7uOSQnd8fLwsy1LdunW1YcMGVa5c2fGct7e3QkJC5OnJV3EAAOBK9GsAANzHJYXu2rVrS5JycnIk5Z6ydujQoTzXi91xxx0lVB4AALhU9GsAANxHsb4yLD4+Xn369NHWrVtls9lkWZak3K8mkXKvFwMAAK5FvwYAwPWKdffyp556SnXq1FFSUpL8/f21fft2rVu3Tm3atFFMTEwJlwgAAIqDfg0AgOsV65Pu2NhYrVmzRpUqVZKHh4c8PT3VoUMHTZo0SU899ZQ2b95c0nUCAIBLRL8GAMD1ivVJt91uV2BgoCSpUqVKOnr0qKTca8h2795dctUBAIBio18DAOB6xfqku1mzZtqyZYvCwsJ0ww03aMqUKfL29tasWbNUt27dkq4RAAAUA/0aAADXK1bofv7553XmzBlJ0ssvv6zbb79dHTt2VMWKFfXpp5+WaIEAAKB46NcAALhesUJ3ZGSk47/Dw8O1a9cunTx5UuXLl3fcERUAALgW/RoAANcrVujOT4UKFUpqUwAAwBD6NQAAV1axbqQGAAAAAAAujtANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQl4budevWqXfv3goNDZXNZtOyZctKZLsxMTFq3bq1fHx8FB4erjlz5uSZc+TIEQ0YMEAVK1aUn5+fmjdvro0bN5bI6wMAUNrQswEAKB6Xhu4zZ86oZcuWeu+990psm/Hx8erVq5e6dOmiuLg4jRw5Uo888oi++uorx5xTp06pffv2KlOmjFatWqUdO3YoKipK5cuXL7E6AAAoTejZAAAUj5crX7xnz57q2bNngc9nZmbqueee0/z583X69Gk1a9ZMkydPVufOnQtcM2PGDIWFhSkqKkqS1LhxY61fv15Tp05VZGSkJGny5MmqWbOmZs+e7VgXFhZWMjsFAEApRM8GAKB43Pqa7uHDhys2NlYLFizQ1q1b1a9fP/Xo0UN79+4tcE1sbKy6d+/uNBYZGanY2FjHz59//rnatGmjfv36KSQkRBEREXr//fcLrSUzM1MpKSlODwAAkIueDQBA/tw2dB86dEizZ8/WokWL1LFjR9WrV09jxoxRhw4dnN7tvtDx48dVpUoVp7EqVaooJSVF6enpkqQDBw5o+vTpql+/vr766is9/vjjeuqpp/Thhx8WuN1JkyYpODjY8ahZs2bJ7CgAAFc5ejYAAAVz6enlhdm2bZvsdrsaNGjgNJ6ZmamKFStKkgICAhzjAwYM0IwZM4q07ZycHLVp00avv/66JCkiIkLbt2/XjBkzNGjQoHzXjB8/XqNGjXL8nJKSQhMHAED0bAAACuO2oTstLU2enp7atGmTPD09nZ4737jj4uIcY0FBQZKkqlWrKjEx0Wl+YmKigoKC5OfnJ0mqVq2amjRp4jSncePGWrx4cYH1+Pj4yMfHp9j7AwBAaUXPBgCgYG4buiMiImS325WUlKSOHTvmOyc8PDzPWLt27bRy5UqnsejoaLVr187xc/v27bV7926nOXv27FHt2rVLoHIAAK4t9GwAAArm0mu609LSFBcX53j3Oz4+XnFxcTp06JAaNGig/v37a+DAgVqyZIni4+O1YcMGTZo0SStWrChwm0OHDtWBAwc0btw47dq1S//617+0cOFCPf300445Tz/9tH744Qe9/vrr2rdvn+bNm6dZs2Zp2LBhpncZAICrEj0bAIDicWno3rhxoyIiIhQRESFJGjVqlCIiIvTCCy9IkmbPnq2BAwdq9OjRatiwoe666y799NNPqlWrVoHbDAsL04oVKxQdHa2WLVsqKipKH3zwgeOrRySpbdu2Wrp0qebPn69mzZrplVde0TvvvKP+/fub3WEAAK5S9GwAAIrHZlmW5eoirkYpKSkKDg5W9xqPy8uD68bcSWrrUFeXgHwkRF58DlzDZtlcXQIukJOeoYTRE5ScnOy4/hnF5+jZVf8uLw9vV5eDvzjdsY6rS0A+lke97eoSUIB2sY+5ugT8hf1shg4MmnTRfu22XxkGAAAAAMDVjtANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgiJerC7jaWWX9ZHn6uLoM/EV2Wd5Lckf+VVJdXQIK8Eu7ua4uARdISc1R+dGurqIU8vSQPDxdXQX+Iq06PdsdjTjcw9UloAC7Onzs6hLwFympOSpfhHn8pgMAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AuCyZmZaeffU31WgVr7Jh+9XutgRFrz1b5PVHjp3T3/5+XBUaHlC5+vt110PHdOBgdr5z/z0vRU07HpR/nf1qeNNBTfv36TxzlqxI032PHVf4Db8qIGy/Gnc4qDETf9PpZHueuaNeOKE2tySoUuMDCgjbr6YdD+qlt35X2pmcItcPAMDV4kr27MQT5zR0bJJqRsTLv85+1W37qx4Zleg0Z+nKNPW474hqtIqXX+19qtU6Xv0eOabtuzKd5sV8f1ae1fYV+Hj9nZOX/odxBXm5ugAAwNVt8MhELV6ephGPllN4WBl9tDBVtw84qtWfVVeHG/wKXZt2Jkfd+h5RckqOxj9VXmW8bHpn1ml1ufuIfo6uqYoVPB1zZ36UrCeeOaG7e5XVyMfKaf2PGRrx/G86m25p3PDyjnlDxyYptKqX+t8TqJrVvbR9Z5bem31aq1af0cava8rP78/3mzfGZarDDb566L5A+fjYFLc9U5OnndbqdemKWVZdHh62kv8DAwDARa5Uz044kq2OdxyRJD32YLBCq3np2PFz2hCX4bTN7buyVL6cp556pJwqVvBUYtI5zV6Qoht7HtZ3y2uoZVMfSVLj+t768N0qeWr65LMURa9N1y2d/S/3j8YoQjcAoNg2bM7Qp8vSNOWFihr9eG7wHdgvUC26JOjZV3/X+i9qFLp++pxk7T2QrR9W1VDbVr6SpB5d/dWiyyG9PeO0Xvu/ipKk9PQcTXjjd93W3V+LPqgmSXp0QLByciy9OvWkHh0QpPLlcpv9wg+qqvNNzs23dQsfDR6RpLlLUvVI/2DH+LrP89ZXr04ZjX3pd23YnKkbr/Mt5p8MAADu5Ur1bEkaOu6EvLykH1c5h/ELTRhVIc/Yw/2DVKv1r5rxYbKmTwmRJFWp7KUBfQPzzH3l7ZOqX7eMox535dLTy9etW6fevXsrNDRUNptNy5YtK5HtxsTEqHXr1vLx8VF4eLjmzJnj9LzdbteECRMUFhYmPz8/1atXT6+88oosyyqR1weAa8Xi5Wny9MwNwOf5+npoyP2Bit2YoYQj+Z9y9tf1bVv5ODXLRvW91bWDnxZ9keYY++b7dP1+KkePPxTstP6JwcE6c9bSiv/9eWrchYFbkvrcFiBJ2rW38HokqXaNMpKU7+noF5r279Nq3umQAsL2q2KjA7o+MkHzlqRedN3ViJ4NAFe3K9Wzd+3N0pdrzmr0E+VVsYKnMjJylJ1d9N/ZIZU85e/nodMphV/qtWFzhvbFZ+uBu/OG8QtlZ1t6OeqkGt6Ue4la5SYHdPMdhy/p1PrL4dLQfebMGbVs2VLvvfdeiW0zPj5evXr1UpcuXRQXF6eRI0fqkUce0VdffeWYM3nyZE2fPl3Tpk3Tzp07NXnyZE2ZMkXvvvtuidUBANeCzdsz1aBuGQUFOreTthG5DTnul6wC1+bkWNq6M0vXtfTJ89z1Eb7a/2u2UtNyG27cttxru9pcMPe6Fr7y8JDitmfm2cZfHU86J0mqVCFv2zt3ztJvv9t19Pg5fR1zVi9M/l2BATZdH1H4u+bvf5KsEc//psYNyujtlytp4pgKatnMRxt+zih03dWKng0AV7cr1bNXf5sbZKtU8tQt/Y6obNgBlQ3br9seOKpfE/IP9qeT7Trxm13bdmbq0dFJSknNUdcOhZ/uPm9x7pvcD9wdUOg8SXrprZN6OeqkOt/kp3++Vkn/91R51azhpc3bCv/3Q0lx6enlPXv2VM+ePQt8PjMzU88995zmz5+v06dPq1mzZpo8ebI6d+5c4JoZM2YoLCxMUVFRkqTGjRtr/fr1mjp1qiIjIyVJ33//ve6880716tVLklSnTh3Nnz9fGzZsKLmdA4BrwPFEu6pVydtKqoXkjh09fq7AtSdP5Sgz03LMdVpf5c/1DcO9dSzJLk9PKaSS81xvb5sqlvfU0cSCX0eSprx3Wp6e0j23523MG7dkqv3thx0/N6xXRss+rKYK5Qs+HU6SVv7vrJo29NbC96sVOq+0oGcDwNXtSvXsvQdyg/XQsUlq08pX82dWUcKRc3o56qRuvfeo4lbXlL+/c/C/qddh7d6fuy6grE3PjSyvhx8IKrAeu93Sws/TdH2Ej8LDvC+y59LK1WfUs5u/Zr4VctG5Jrj13cuHDx+u2NhYLViwQFu3blW/fv3Uo0cP7d27t8A1sbGx6t69u9NYZGSkYmNjHT/fdNNNWr16tfbs2SNJ2rJli9avX3/Rf0ykpKQ4PQDgWpeekSNv77w3G/P1zR3LyCj4dLL0jNx3xH188q4/P5b+x/r0jBx5l8n/pma+Pjalpxf8OvOWpOo/81I0amg51a+btzE3aeCtrz4N1ZLZVTV2WDmV9fdQ2pmLnwZXLthDh4+d009xpfOT7UtFzwYA93alenba2dy5VUO8tPyTarr3jkCNfry8Zr4Vov2/Zmve0ryXYf37nRCtnFdN096orMb1vZWeYcleyFVeq79NV+IJe5FOLZekckEe2rE7S3sPFPxpvklueyO1Q4cOafbs2Tp06JBCQ0MlSWPGjNGXX36p2bNn6/XXX8933fHjx1WlivOd7apUqaKUlBSlp6fLz89Pzz77rFJSUtSoUSN5enrKbrfrtddeU//+/QusZ9KkSXrppZdKbgcBoBTw8/VQVlbeJn2+cZ9v5AWtlXK/vuRC58f8/ljv5+uhrAKuB8vItOTnl//rfPtDuh4dnaRbO/vr1Wcr5jsnKNBD3W/OvQ78zh4BmrckVX0eOqaNX9d03DU1P+OGldfqb8/qxp6HFR5WRrd08tf9fQLU/vrCT4crjejZAOD+rmTPlqR+dwQ4fQtIv94BGvRkomI3Zjjd1FSS2rX5s3fed2eAmt58SJL05ouV8q1n3pJUeXpK99558VPLJWni2IrqM/iYGrU/pGaNvBXZxV8D+gaqRZOC+3xJcttPurdt2ya73a4GDRooICDA8Vi7dq32798vSU7jQ4cOLfK2Fy5cqLlz52revHn6+eef9eGHH+qtt97Shx9+WOCa8ePHKzk52fFISEi47H0EgKtd1SqeOpbPqd3H/riGOrRqwe/tVijvIR8fm2Ou0/pE5/XVQjxlt0tJvznPzcqy9Pspu0LzOV1uyy+ZuuuhY2rW0FuLPqgqL6+iff3X3beVlSR9uqzwG6I1buCtnetra96MKmp/va+WrEjTzXce0cQ3fy/S65Qm9GwAcH9XqmeHVsm9PCukkvNlWp6euZeEnTpd+A3SypfzVJcOfgXemDQ9PUfLVqWpW0d/ValctM+Qb27np72xtfXB1BA1beStf89LUZtbE/TB3OQirb9cbvtJd1pamjw9PbVp0yZ5ejofsICA3Hc04uLiHGNBQbnn/FetWlWJic5fup6YmKigoCD5+eW+gzJ27Fg9++yzuu+++yRJzZs318GDBzVp0iQNGjQo33p8fHzk43Nl3gkBgKtFq6Y+ivkuXSmpOU43Zjl/M7FWTQu+zsrDw6bmjby1aUvem5j8+HOG6tb2UmBA7jZbNsv9/btxS6Zu6/Zn69q4JUM5OVLLC15n/6/Zuu2Bowqp5Knln4QqoGzR32POzLKUkyMlpxb+jwJJKuvvob/dGai/3RmorCxL9zx8TK//45SefbK8fH3d9n3tEkfPBgD3d6V6dusWub9/L7xGPCvL0m8n7apcsfB7pkhSerql5ALuXv7512eUmmYV6QZqf1WhvKcG3xekwfcFKe1MjjrfdUQvR53M86m7CW77L4KIiAjZ7XYlJSUpPDzc6VG1alVJchoLCcm9KL5du3ZavXq107aio6PVrl07x89nz56Vh4fzrnt6eion5+L/wAIA/Ome2wNkt+feyfu8zExLcz5N1Q2tfVSzehnH+KHD2dq1NyvP+p/iMrXxL9dF796XpW++S1ff3n82067t/VShvIdmfOj8jvSMD1Pk72dTr+5lHWPHk86px31H5OEhrZofqsqV8m/up5Pt+X6Fyb/n5V7/e13Lwu9e/vtJ54vNvL1tatLAW5YlZRd+X7dSh54NAO7vSvXszjf5K6SSp+YtSVVGxp+/q+d8miK7Xere6c+v9rzwDDZJ+jUhW2vWn83zjSXnzV+aJn8/m+PrQIviwp4dUNZD9cLK5Hu6vAku/aQ7LS1N+/btc/wcHx+vuLg4VahQQQ0aNFD//v01cOBARUVFKSIiQidOnNDq1avVokULx11MLzR06FBNmzZN48aN05AhQ7RmzRotXLhQK1ascMzp3bu3XnvtNdWqVUtNmzbV5s2b9fbbb2vIkCHG9xkASpMbWvuqb+8A/d/rvyvpN7vq1Smjjxel6teEbL0fVd1p7kNPJWptbIbsx8IdY48/FKwP5qao94PHNPrxcirjZdPUmadVpbKnRj1W3jHPz89DL4+rqOHjT+jeR4/p1s7+Wv9jhuYuTtWrz1ZwutP4bQ8c1YGD5zR2WDmt/zFD63/88x8HVSp76pY/mn3M9+ka+fxvuuf2sgqv663sLEvf/piupSvPqE1LHw24p/Cbs/S476iqhHiqfVtfhVT20q69WXpvdrJu6+7veLe/NKFnA8DV7Ur1bB8fmyZPqKjBI5LUuc8RDegbqENHzumfH5xWxxt8HZdxSVLLLgnq2tFPrZr6qFywh/bFZ+s/81OUfU56/bm892I5ecquL9ec0d29Ai7pLLZmnQ6p001+uq6Fj8qX89CmLZlavDxNw4aY/5RbkmyWZV2ZeJ+PmJgYdenSJc/4oEGDNGfOHGVnZ+vVV1/VRx99pCNHjqhSpUq68cYb9dJLL6l58+aFbvfpp5/Wjh07VKNGDU2YMEEPPfSQ4/nU1FRNmDBBS5cuVVJSkkJDQ3X//ffrhRdekLf3xW85L0kpKSkKDg5Wt4aj5OXJKWzu5FRE/jdLgmul3FP49bFwnV/azb2s9RkZOXphyknNXZyqU8k5atHYWy+Nq6DILmWd5nW9+3CeBi5Jh4+e06gXTyh6bbpycix1uslPb79UKd+vAHn/k2RNnXFa8QnZqhlaRsMGB+upR4Nls/15vbZntX151p3XqZ2v1iypISn3FPRX3j6p7zak61iiXZakerXL6J7by2rME+VV1r/wZj7r42TNX5KqX3ZnKe2spRrVvNTntrJ6bmSFPN+BeqlSUnNUvsEBJScnO07FdrXS0LO7Vx8qLw96tjs53K+2q0tAPlreu93VJaAAH9Ved1nrr2TPXrAsVVOmndKufdkqF+Shvr0D9Nr4ik5vTL/01u9a+b+z2n8w93u+Qyp5quONfhr/VHk1b5z39/XMj5L1xDMntOzDaup9a9k8zxfk9XdO6ouvz2jPgWxlZlqqXcNLA/oGaswT5VWmgG9HKYqi9muXhu6rGaHbfRG63ROh231dbuhGyXPH0H01I3S7L0K3eyJ0u6/LDd0oWUXt16Xv/DcAAAAAANwEoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QAAAAAAGELoBgAAAADAEEI3AAAAAACGELoBAAAAADCE0A0AAAAAgCGEbgAAAAAADCF0AwAAAABgCKEbAAAAAABDCN0AAAAAABhC6AYAAAAAwBBCNwAAAAAAhhC6AQAAAAAwhNANAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAY4uXqAq5WlmVJks7ZM11cCS5kz85wdQnIh/0sf1fcVUpqjqtLwAVS0nKPyfleg8vj6Nk5WS6uBBeyZ9Kz3VH2Gf6uuCt6tnspar+2WXT0Yjl8+LBq1qzp6jIAAKVYQkKCatSo4eoyrnr0bACASRfr14TuYsrJydHRo0cVGBgom83m6nIuS0pKimrWrKmEhAQFBQW5uhz8gePinjgu7qm0HRfLspSamqrQ0FB5eHAl2OWiZ8M0jot74ri4p9J0XIrarzm9vJg8PDxK3acPQUFBV/3/+KURx8U9cVzcU2k6LsHBwa4uodSgZ+NK4bi4J46Leyotx6Uo/Zq3zwEAAAAAMITQDQAAAACAIYRuyMfHRy+++KJ8fHxcXQr+guPinjgu7onjgmsF/6+7J46Le+K4uKdr8bhwIzUAAAAAAAzhk24AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQui+yq1bt069e/dWaGiobDabli1bViLbjYmJUevWreXj46Pw8HDNmTPH6flJkyapbdu2CgwMVEhIiO666y7t3r27RF67NHDVcfmrN954QzabTSNHjiyR1y4NXHlcjhw5ogEDBqhixYry8/NT8+bNtXHjxhJ5/audq46L3W7XhAkTFBYWJj8/P9WrV0+vvPKKuL8oTKBfuyf6tXuiX7svenbxELqvcmfOnFHLli313nvvldg24+Pj1atXL3Xp0kVxcXEaOXKkHnnkEX311VeOOWvXrtWwYcP0ww8/KDo6WtnZ2br11lt15syZEqvjauaq43LeTz/9pJkzZ6pFixYl9vqlgauOy6lTp9S+fXuVKVNGq1at0o4dOxQVFaXy5cuXWB1XM1cdl8mTJ2v69OmaNm2adu7cqcmTJ2vKlCl69913S6wO4Dz6tXuiX7sn+rX7omcXk4VSQ5K1dOlSp7GMjAxr9OjRVmhoqOXv729df/311jfffFPodsaNG2c1bdrUaexvf/ubFRkZWeCapKQkS5K1du3a4pZfal3p45KammrVr1/fio6Otjp16mSNGDGiBPai9LmSx+WZZ56xOnToUFKll2pX8rj06tXLGjJkiNOcu+++2+rfv/9l7QNwMfRr90S/dk/0a/dFzy46Puku5YYPH67Y2FgtWLBAW7duVb9+/dSjRw/t3bu3wDWxsbHq3r2701hkZKRiY2MLXJOcnCxJqlChQskUXsqZPC7Dhg1Tr1698szFxZk6Lp9//rnatGmjfv36KSQkRBEREXr//feN7UdpY+q43HTTTVq9erX27NkjSdqyZYvWr1+vnj17mtkRoBD0a/dEv3ZP9Gv3Rc8ugKtTP0qOLni36eDBg5anp6d15MgRp3ndunWzxo8fX+B26tevb73++utOYytWrLAkWWfPns0z3263W7169bLat29/eTtQSl3J4zJ//nyrWbNmVnp6umVZFu+cF+JKHhcfHx/Lx8fHGj9+vPXzzz9bM2fOtHx9fa05c+aU3A6VElfyuNjtduuZZ56xbDab5eXlZdlstjxrABPo1+6Jfu2e6Nfui55ddF5XPubjStm2bZvsdrsaNGjgNJ6ZmamKFStKkgICAhzjAwYM0IwZMy75dYYNG6bt27dr/fr1l1fwNcLUcUlISNCIESMUHR0tX1/fki36GmDy70tOTo7atGmj119/XZIUERGh7du3a8aMGRo0aFAJ7UHpZPK4LFy4UHPnztW8efPUtGlTx3VkoaGhHBdcUfRr90S/dk/0a/dFzy4YobsUS0tLk6enpzZt2iRPT0+n587/Dx8XF+cYCwoKkiRVrVpViYmJTvMTExMVFBQkPz8/p/Hhw4dr+fLlWrdunWrUqGFgL0ofU8dl06ZNSkpKUuvWrR3P2+12rVu3TtOmTVNmZmae18OfTP59qVatmpo0aeI0p3Hjxlq8eHFJ70apY/K4jB07Vs8++6zuu+8+SVLz5s118OBBTZo06apo4Cg96NfuiX7tnujX7oueXTBCdykWEREhu92upKQkdezYMd854eHhecbatWunlStXOo1FR0erXbt2jp8ty9KTTz6ppUuXKiYmRmFhYSVbfClm6rh069ZN27Ztc3p+8ODBatSokZ555hka+EWY/PvSvn37PF/Rs2fPHtWuXbsEKi/dTB6Xs2fPysPD+dYmnp6eysnJKYHKgaKjX7sn+rV7ol+7L3p2IVx9fjsuT2pqqrV582Zr8+bNliTr7bfftjZv3mwdPHjQsizL6t+/v1WnTh1r8eLF1oEDB6wff/zRev31163ly5cXuM0DBw5Y/v7+1tixY62dO3da7733nuXp6Wl9+eWXjjmPP/64FRwcbMXExFjHjh1zPPK7huxa5KrjciGuEXPmquOyYcMGy8vLy3rttdesvXv3WnPnzrX8/f2tTz75xPg+Xw1cdVwGDRpkVa9e3Vq+fLkVHx9vLVmyxKpUqZI1btw44/uMaw/92j3Rr90T/dp90bOLh9B9lfvmm28sSXkegwYNsizLsrKysqwXXnjBqlOnjlWmTBmrWrVqVp8+faytW7dedLutWrWyvL29rbp161qzZ892ej6/15SUZ961ylXH5UI0cWeuPC5ffPGF1axZM8vHx8dq1KiRNWvWLAN7eHVy1XFJSUmxRowYYdWqVcvy9fW16tataz333HNWZmamoT3FtYx+7Z7o1+6Jfu2+6NnFY7MsyyqZz8wBAAAAAMBf8T3dAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGAIoRsAAAAAAEMI3QDchs1m07Jly1xdBgAAKAT9Grg0hG4AJY5mDACA+6NfA1cGoRtAicrKynJ1CQAA4CLo18CVQ+gGcFk6d+6s4cOHa+TIkapUqZJ8fHwkSX369JHNZlOdOnUcc6dPn6569erJ29tbDRs21Mcff+yiqgEAuLbQrwHXIXQDuGwffvihvL299d133+mHH36QJM2ePVvHjh3TTz/9JElaunSpRowYodGjR2v79u167LHHNHjwYH3zzTeuLB0AgGsG/RpwDZtlWZariwBw9ercubNSUlL0888/O8ZsNpuWLl2qu+66yzHWvn17NW3aVLNmzXKM3XvvvTpz5oxWrFhR4DoAAHD56NeA6/BJN4DLdt111110zs6dO9W+fXunsfbt22vnzp2mygIAAH9BvwZcg9AN4LKVLVvW1SUAAICLoF8DrkHoBlDiypQpI7vd7jTWuHFjfffdd05j3333nZo0aXIlSwMAAH+gXwNXhperCwBQ+tSpU0erV69W+/bt5ePjo/Lly2vs2LG69957FRERoe7du+uLL77QkiVL9L///c/V5QIAcE2iXwNXBp90AyhxUVFRio6OVs2aNRURESFJuuuuu/SPf/xDb731lpo2baqZM2dq9uzZ6ty5s2uLBQDgGkW/Bq4M7l4OAAAAAIAhfNINAAAAAIAhhG4AAAAAAAwhdAMAAAAAYAihGwAAAAAAQwjdAAAAAAAYQugGAAAAAMAQQjcAAAAAAIYQugEAAAAAMITQDQAAAACAIYRuAAAAAAAMIXQDAAAAAGDI/wMc+2TI/0ZZsgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# time how long it takes to solve the model with different tolerances\n",
    "models = [pybamm.lithium_ion.SPM(), pybamm.lithium_ion.DFN()]\n",
    "model_names = [\"SPM\", \"DFN\"]\n",
    "atols = [1e-2, 1e-4, 1e-6, 1e-8]\n",
    "rtols = [1e-2, 1e-4, 1e-6, 1e-8]\n",
    "results = np.zeros((len(models), len(atols), len(rtols)))\n",
    "for imodel, model in enumerate(models):\n",
    "    for iatol, atol in enumerate(atols):\n",
    "        for irtol, rtol in enumerate(rtols):\n",
    "            solver = pybamm.IDAKLUSolver(atol=atol, rtol=rtol)\n",
    "            sim = pybamm.Simulation(model, solver=solver)\n",
    "            sol = sim.solve([0, 3600])\n",
    "            start_time = time.perf_counter()\n",
    "            sol = sim.solve([0, 3600])\n",
    "            end_time = time.perf_counter()\n",
    "            results[imodel, iatol, irtol] = end_time - start_time\n",
    "\n",
    "# plot results in a separate plot for each model\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n",
    "for imodel, _model in enumerate(models):\n",
    "    ax[imodel].imshow(\n",
    "        results[imodel],\n",
    "        vmin=results[imodel].min(),\n",
    "        vmax=results[imodel].max(),\n",
    "        cmap=\"viridis\",\n",
    "    )\n",
    "    ax[imodel].set_xticks(range(len(rtols)))\n",
    "    ax[imodel].set_xticklabels([f\"{rtol:.0e}\" for rtol in rtols])\n",
    "    ax[imodel].set_xlabel(\"rtol\")\n",
    "    ax[imodel].set_yticks(range(len(atols)))\n",
    "    ax[imodel].set_yticklabels([f\"{atol:.0e}\" for atol in atols])\n",
    "    ax[imodel].set_ylabel(\"atol\")\n",
    "    ax[imodel].set_title(model_names[imodel])\n",
    "    ax[imodel].text(\n",
    "        0,\n",
    "        0,\n",
    "        f\"{results[imodel].min():.4f} s\",\n",
    "        ha=\"center\",\n",
    "        va=\"center\",\n",
    "        color=\"white\",\n",
    "        fontsize=12,\n",
    "    )\n",
    "    ax[imodel].text(\n",
    "        3,\n",
    "        3,\n",
    "        f\"{results[imodel].max():.4f} s\",\n",
    "        ha=\"center\",\n",
    "        va=\"center\",\n",
    "        color=\"black\",\n",
    "        fontsize=12,\n",
    "    )\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, the choice of tolerances is a trade-off between accuracy and computational time. It can often be difficult to know what tolerances to use, but a good starting point is to use the default tolerances and then adjust them if the solver is too slow or too inaccurate. \n",
    "\n",
    "You can also use a convergence study to determine the best tolerances for your problem, using either an analytical solution (if available) or by comparing the solution to a high-tolerance solution. For example, lets compare the solution of a high-tolerance SPM model to a lower-tolerance SPM model while varying the tolerances. We'll pick some goal relative and absolute tolerances of $10^{-4}$ and $10^{-6}$ respectively and use a similar error norm as the solve. We'll then then vary the tolerances from $10^{-2}$ to $10^{-6}$ and see when the solution is within the set goal tolerances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG1CAYAAAAV2Js8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABQRklEQVR4nO3de1xUdeL/8deAXLwAiiiIgnhBExUoBbIyNdlVK91yK7e0UNPdyswi29Vfba3frWyrLXNl1y0tu7ibWZt2dSvUzLzgJbxfQFHxAogoV7nNnN8fJhspwijDGZj38/HgsTtnzsy855Myb8/5fM5YDMMwEBEREXFBbmYHEBERETGLipCIiIi4LBUhERERcVkqQiIiIuKyVIRERETEZakIiYiIiMtSERIRERGXpSIkIiIiLquZ2QGcnc1m4/jx4/j4+GCxWMyOIyIiInVgGAaFhYUEBwfj5lbzcR8VoVocP36ckJAQs2OIiIjIZcjMzKRTp0413q8iVAsfHx/g3ED6+vqanEZERETqoqCggJCQkKrP8ZqoCNXi/OkwX19fFSEREZFGprZpLZosLSIiIi5LRUhERERcloqQiIiIuCwVIREREXFZKkIiIiLislSERERExGWpCImIiIjLUhESERERl+USRej222+nTZs23HHHHWZHEREREcBqM1h/4BTLU4+x/sAprDbDlBwucWXpadOmMXHiRN5++22zo4iIiLi8FTtPMOvT3ZzIL63a1sHPm2dGRjC8T4cGzeISR4QGDx5c63eNiIiIiOOt2HmCB9/bWq0EAWTll/Lge1tZsfNEg+Zx+iK0Zs0aRo4cSXBwMBaLhWXLll2wT1JSEmFhYXh7exMXF0dKSkrDBxUREZFLstoMZn26m4udBDu/bdanuxv0NJnTF6Hi4mKioqJISkq66P1LliwhMTGRZ555hq1btxIVFcWwYcPIyclp4KQiIiJyKSkZeRccCfopAziRX0pKRl6DZXL6OUIjRoxgxIgRNd7/yiuvMHnyZCZMmADA/Pnz+fzzz3nzzTeZMWOG3a9XVlZGWVlZ1e2CggL7Q4uIiMgFcgprLkGXs199cPojQpdSXl7Oli1biI+Pr9rm5uZGfHw869evv6znnD17Nn5+flU/ISEh9RVXRETEZeWfrWDN/pN12re9j7eD0/xPoy5Cubm5WK1WAgMDq20PDAwkKyur6nZ8fDx33nknX3zxBZ06dbpkSZo5cyb5+flVP5mZmQ7LLyIi0tTll1Twytf7ueEvK/lo67FL7mvh3Oqx2C7+DROORnBqrD588803dd7Xy8sLLy8vkpKSSEpKwmq1OjCZiIhI03SmpJyFazNY9P0hCssqAegR2Iobe7Rj4XcZANUmTVt+/N9nRkbg7mahoTTqIhQQEIC7uzvZ2dnVtmdnZxMUFHRFzz1lyhSmTJlCQUEBfn5+V/RcIiIiruJ08Y8FaN0hin4sQD0DfZgWH87w3kG4uVno37nNBdcRCjLpOkKNugh5enrSr18/kpOTue222wCw2WwkJyfz8MMPmxtORETEheQVl7Pgu4O8ve4QxeXnzqZcFeTDtKHhDPuxAJ03vE8HfhERREpGHjmFpbT3OXc6rCGPBJ3n9EWoqKiI9PT0qtsZGRmkpqbi7+9PaGgoiYmJJCQk0L9/f2JjY5kzZw7FxcVVq8gul06NiYiI1O5UURlvfJfBO+sPUfJjAYro4MsjQ8P5ZURgtQL0U+5uFgZ0a9uQUS/KYhiGOV/uUUerV69myJAhF2xPSEhg0aJFAMybN4+XXnqJrKwsoqOjmTt3LnFxcfXy+udPjeXn5+Pr61svzykiItLY5RaV8caag7y74XBVAeod7Mu0oeH8IiIQi6Xhj+78VF0/v52+CJlNRUhEROR/ThaW8fqaA7y34QhnK84VoD4dfXl0aA+G9mpvegE6r66f305/akxERETMl1NYyuvfHuS9jYcprbABENnJj2lDw7npKucpQPZSEaqB5giJiIhATkEp8789yOKNhymrPFeAokJa8+jQcAb3bNdoC9B5OjVWC50aExERV5RdUMo/Vh/g3ylHqgrQ1aGtmTY0nEE9nL8A6dSYiIiI2C0rv5R/rE7n35syKf+xAPXr3IZpQ8MZGB7g9AXIXipCIiIiwon8s/xj9QHeT8mk3HquAPXv3IZH43twffe2Ta4AnaciVAPNERIREVdw/MxZ/r46nQ82Ha0qQLFh/kyLD+e6bk23AJ2nOUK10BwhERFpio6eLuHvqw+wdHMmFdZzVSCuy7kCNKBr4y9AmiMkIiIiF8jMO1eAPtzyvwJ0bVd/pg3t4RRXem5oKkIiIiIuIDOvhKRV6Xy45SiVtnMF6LpubZk2NJy4rq5XgM5TEaqB5giJiEhTcORUCfNWpfGfrceqCtAN3QOYFh9OTJi/yenMpzlCtdAcIRERaYwO5RYzb1U6H/9wDOuPBWhgeACPxofTr3PTL0CaIyQiIuKCMnKLmbcynWWp/ytAN/Zox7Sh4fTr3MbkdM5HRUhERKQJOHiyqKoA/dh/GNzzXAG6OlQFqCYqQiIiIo1Yek4R81am8cm241UF6Kar2vPI0HCiQ1qbmq0xUBESERFphNJzCpmbnM6n249zfrZvfK9zBSiyU2tTszUmKkI10KoxERFxRvuzC5mbnMbnO078pAAFMm1oOH07+ZkbrhHSqrFaaNWYiIg4g31ZhcxdmcYXPylAv4wI5JGh4fTpqAL0c1o1JiIi0gTszSpgbnIaX+zIqto2vHcQjwwNJyJY/0C/UipCIiIiTmj38XMFaMWu/xWgm/sGMfWmcHp1UAGqLypCIiIiTmTnsXzmJqfx1e5sACwWuLlvBx65KZyeQT4mp2t6VIREREScwM5j+byWnMbXPylAt/TtwCNDw+kRqALkKCpCIiIiJtpxNJ/XkvfzzZ4c4FwBGhkZzNSbuhOuAuRwKkI10PJ5ERFxpG2ZZ3gtOY2Ve88VIDcLjIoK5uGbwunevpXJ6VyHls/XQsvnRUSkPv1w5DSvJaexet9J4FwBui26I1Nu6k63dipA9UXL50VERJzI1iOnee2bNL7df64AubtZ+FV0MA8P6U5XFSDTqAiJiIg40JbDecz5Jo3v0nKBcwXo9qs78vCQ7oQFtDQ5nagIiYiIOMCmQ3m89k0aa9P/V4B+fU1HpgzpTue2KkDOQkVIRESkHm08eIrXktNYd+AUAM3cLNzRrxNThnQnxL+Fyenk51SERERE6sGGg6d47Zs01h/8XwG6s38nHhqsAuTMVIREREQuk2EYrP+xAG3MyAPAw93Cnf1DeGhwNzq1UQFydipCIiIidjIMg/UHTjHnmzRSDp0rQJ7ubtwV04kHB3enY+vmJieUulIREhERqSPDMPg+/RSvJe9n06HTwLkC9JvYEB4Y1I1gFaBGR0WoBrqytIiInGcYBt+l5fJachpbDv9YgJq5cXdMCA8M7kYHPxWgxkpXlq6FriwtIuK6DMPg2/0neS05jR+OnAHAq5kbd8eG8uDgbgT6epsbUGqkK0uLiIhcJsMwWL3vXAFKzTwDnCtAY+M688CgrrRXAWoyVIRERER+ZBgGq/bl8No3aWw7mg+At4cb4+I689tBXWnvowLU1KgIiYiIyzMMg+Q9Ocxdmcb2nxSge6/tzG9v7EY7Hy+TE4qjqAiJiIjLMgyDr3dnM3dlGjuPFQDQ3MOd+wZ0ZvKNXQlopQLU1KkIiYiIy7HZDL7anc3c5DR2nzhXgFp4unPfgDAmD+xCWxUgl6EiJCIiLuNcAcpizjdp7M0qBKClpzsJ14UxaWBX/Ft6mpxQGpqKkIiINHk2m8GKXVnMTf5fAWrl1YyE6zoz6YautFEBclkqQiIi0mTZbAZf7DzB35LT2Zd9rgD5eDVj/PVh3H9DF1q3UAFydSpCIiLS5FhtBp/vOMHfktNIyykCzhWgCTd04f7ru+DXwsPkhOIsVIRERKTJsNoMPtt+nL+tTCf9fAHybsbE67swUQVILsIlitBnn33G448/js1m4w9/+AOTJk0yO5KIiNQjq83g023H+dvKNA6cLAbA17sZ99/QlfHXh+HXXAVILq7JF6HKykoSExNZtWoVfn5+9OvXj9tvv522bduaHU1ERK5QpdXGJ9uOM29lOgdzzxUgv+YeTLqhCwnXh+HrrQIkl9bki1BKSgq9e/emY8eOAIwYMYKvvvqKu+++2+RkIiJyuSqtNpalHidpVToZPxag1i08mDywK/cN6IyPCpDUkZvZAWqzZs0aRo4cSXBwMBaLhWXLll2wT1JSEmFhYXh7exMXF0dKSkrVfcePH68qQQAdO3bk2LFjDRFdRETqWaXVxtLNmQx95VumL91GRm4xbVp48PvhPVn7h5uYMqS7SpDYxemPCBUXFxMVFcXEiRMZPXr0BfcvWbKExMRE5s+fT1xcHHPmzGHYsGHs27eP9u3bm5BYRETqW4XVxsdbjzFvVTpH8koA8G/pWXUEqKWX03+ciZNy+j85I0aMYMSIETXe/8orrzB58mQmTJgAwPz58/n888958803mTFjBsHBwdWOAB07dozY2Ngan6+srIyysrKq2wUFBfXwLkRE5HJUWG18tOUoSavTycw7C0Dblp789saujLtWBUiuXKP+E1ReXs6WLVuYOXNm1TY3Nzfi4+NZv349ALGxsezcuZNjx47h5+fHl19+yR//+Mcan3P27NnMmjXL4dlFRKRm5ZU2PtxylKRV6Rw7c64ABbTy5Hc3dmPstaG08GzUH1/iRBr1n6Tc3FysViuBgYHVtgcGBrJ3714AmjVrxl//+leGDBmCzWbj97///SVXjM2cOZPExMSq2wUFBYSEhDjmDYiISDVllVaWbj7KP1Yf+EkB8uKBQV0ZG9eZ5p7uJieUpqZRF6G6GjVqFKNGjarTvl5eXnh56VuHRUQaUlmllQ82H+Ufq9I5nl8KQDsfLx4Y1I17YkNVgMRhGnURCggIwN3dnezs7Grbs7OzCQoKuqLnTkpKIikpCavVekXPIyIiNSutsPLB5kz+sfoAJ34sQO19vHhwcDfujg3F20MFSByrURchT09P+vXrR3JyMrfddhsANpuN5ORkHn744St67ilTpjBlyhQKCgrw8/Orh7QiInJeaYWV91OO8I9vD5BdcG6BSpCvNw8O7saYmBAVIGkwTl+EioqKSE9Pr7qdkZFBamoq/v7+hIaGkpiYSEJCAv379yc2NpY5c+ZQXFxctYpMREScR2mFlX9tPML8bw+QU3iuAHXw8+ahwd24s78KkDQ8py9CmzdvZsiQIVW3z09kTkhIYNGiRYwZM4aTJ0/y9NNPk5WVRXR0NCtWrLhgArW9dGpMRKT+lFZYWfxjATr5YwEK9vPmwSHduat/J7yaqQCJOSyGYRhmh3Bm50+N5efn4+vra3YcEZFG5Wy5lcUbDzP/24PkFp0rQB1bN+ehId24o58KkDhOXT+/nf6IkIiIND4l5ZW8t+Ewr685SG5ROQCd2jRnypDu/PqaTng2c/pveBIXoSJUA50aExGxX3FZJe9uOMwbaw5yqvh/BWjqTd0ZfU0nPNxVgMS56NRYLXRqTESkdsVllbyz/jBvfHeQvB8LUKh/Cx4e0p3br+moAiQNTqfGRESkXlhtBikZeeQUltLex5vYLv64u1kAKCqr5O11h1jw3UFOl1QA0LntuQJ029UqQOL8VIRERKRGK3aeYNanu6sudgjnlrv/flhPjp05y4K1GZz5sQB1CWjJw0O686voYJqpAEkjoSJUA80REhFXt2LnCR58bys/nz9xIr+Uxz7YVnW7a0BLpg7tzshIFSBpfDRHqBaaIyQirshqM7jhLyurHQn6OXc3Cy/fEcmo6I5Vp8pEnEVdP79V3UVE5AIpGXmXLEFwriwF+TVXCZJGTUVIREQukFN46RJk734izkpFSERELtDex7te9xNxVipCNUhKSiIiIoKYmBizo4iINLhQ/xa4W2o+5WXh3Oqx2C7+DRdKxAE0WboWmiwtIq4mt6iMMf9cz4GTxRe9/3w9+se4axjep0PDBROxgyZLi4iI3c6UlDNuwUYOnCwm2M+bZ2/rQwe/6qe/gvy8VYKkydB1hEREBICC0gruezOFvVmFtPPxYvHka+kS0JK7Y0NrvLK0SGOnIiQiIpSUVzLxrU1sP5pPmxYeLJ4UR5eAlsC56wUN6NbW5IQijqFTYzXQZGkRcRWlFVYmvb2ZzYdP4+PdjHfvj6NHoI/ZsUQahCZL10KTpUWkKSuvtPG7dzezat9JWnq68+6kOK4JbWN2LJErpsnSIiJySZVWG4/8+wdW7TuJt4cbC8fHqASJy1EREhFxQVabwfSl21ixKwtPdzdev7c/13bVPCBxPSpCIiIuxjAMnvx4B8tSj9PMzcLfx17DjT3amR1LxBQqQiIiLsQwDGZ9upv3N2XiZoE5v4kmPiLQ7FgiplEREhFxEYZh8JcV+1i07hAAL90Rxa2RweaGEjGZilANtHxeRJqav61MZ/63BwB47vY+/LpfJ5MTiZhPy+droeXzItIUvL7mAM9/sReAP94awf03dDE5kYhjafm8iIgA8O76Q1Ul6IlhPVWCRH5CRUhEpAn7YHMmf1y+C4CHh3RnypDuJicScS4qQiIiTdTy1GP84aPtANx/Qxce/2UPkxOJOB8VIRGRJmjFziwSP9iGYcDYuFCeuqUXFou+MV7k51SERESamNX7cpj6761YbQajr+nIn3/VRyVIpAYqQiIiTci6A7n87t0tVFgNbonswIu/jsTNTSVIpCYqQiIiTcSWw3lMenszZZU24nsFMmdMNM3c9Wte5FL0N0REpAnYcTSf8W9uoqTcysDwAObdczUeKkEitdLfEhGRRm5vVgH3vrmRwrJKYrv48/q9/fH2cDc7lkijoCJUA33Fhog0BgdOFjFuwUbOlFQQHdKaN8fH0NxTJUikrvQVG7XQV2yIiLM6cqqEu/65nqyCUnoH+/Kvydfi19zD7FgiTkFfsSEi0oQdP3OWexZsIKuglPD2rXj3/jiVIJHLoCIkItLI5BSWMnbBRo6ePkuXgJYsnhSHf0tPs2OJNEoqQiIijUhecTnjFmwkI7eYjq2bs3hSHO19vc2OJdJoqQiJiDQS+WcruHfhRvZnFxHo68W/J19LcOvmZscSadRUhEREGoGiskrGv5XCruMFtG3pyeJJ1xLatoXZsUQaPRUhEREnd7bcyv2LNvHDkTP4NffgvUlxdG/fyuxYIk2CipCIiBMrq7Ty23c3szEjDx+vZrx7fyy9OuhSHiL1RUVIRMRJVVhtTFn8A9+l5dLcw523JsQQ2am12bFEmhQVIRERJ2S1GTy6JJVv9mTj2cyNhQn96R/mb3YskSZHRUhExMnYbAa//3A7n28/gYe7hX+O68d13QPMjiXSJLlEEbr99ttp06YNd9xxh9lRREQuyTAM/rh8Jx9tPYq7m4W/3X01Q65qb3YskSbLJYrQtGnTeOedd8yOISJySYZh8Ozne1i88QgWC7xyVxTD+3QwO5ZIk+YSRWjw4MH4+PiYHUNE5JJe+Xo/C9dmAPCX0ZH8KrqjyYlEmj7Ti9CaNWsYOXIkwcHBWCwWli1bdsE+SUlJhIWF4e3tTVxcHCkpKQ0fVETEgZJWpfO3lekAzBrVm7tiQkxOJOIampkdoLi4mKioKCZOnMjo0aMvuH/JkiUkJiYyf/584uLimDNnDsOGDWPfvn20b3/uvHl0dDSVlZUXPParr74iODjYrjxlZWWUlZVV3S4oKLDzHYmI2OfNtRm89N99AMwccRUJ14WZG0jEhZhehEaMGMGIESNqvP+VV15h8uTJTJgwAYD58+fz+eef8+abbzJjxgwAUlNT6y3P7NmzmTVrVr09n4jIpfxr4xH+77PdADwaH87vBnUzOZGIazH91NillJeXs2XLFuLj46u2ubm5ER8fz/r16x3ymjNnziQ/P7/qJzMz0yGvIyLyn61HeXLZDgB+N6gr04aGm5xIxPWYfkToUnJzc7FarQQGBlbbHhgYyN69e+v8PPHx8Wzbto3i4mI6derE0qVLGTBgwEX39fLywsvL64pyi4jU5vPtJ5i+dBuGAQkDOjNj+FVYLBazY4m4HKcuQvXlm2++sfsxSUlJJCUlYbVaHZBIRFxZ8p5spr3/AzYDxvQP4ZmRvVWCREzi1KfGAgICcHd3Jzs7u9r27OxsgoKCHPraU6ZMYffu3WzatMmhryMiruW7tJM8+N5WKm0Gv4oO5vnRfXFzUwkSMYtTFyFPT0/69etHcnJy1TabzUZycnKNp7ZERJzVxoOnmPzOZsqtNob3DuKvd0bhrhIkYirTT40VFRWRnp5edTsjI4PU1FT8/f0JDQ0lMTGRhIQE+vfvT2xsLHPmzKG4uLhqFZmj6NSYiNSnH46cZuKiTZRW2BjSsx1z776aZu5O/W9REZdgMQzDMDPA6tWrGTJkyAXbExISWLRoEQDz5s3jpZdeIisri+joaObOnUtcXFyD5CsoKMDPz4/8/Hx8fX0b5DVFpGnZdTyfu1/fQEFpJdd1a8ub42Pw9nA3O5ZIk1bXz2/Ti5CzUxESkSuRll3ImNc3kFdcTv/ObXh7YiwtvUw/GC/S5NX181vHZWuQlJREREQEMTExZkcRkUYqI7eYexZsJK+4nMhOfrw5IUYlSMTJ6IhQLXRESEQux9HTJdw1fz3H80u5KsiH9397La1beJodS8Rl6IiQiIhJsvJLueeNjRzPL6Vbu5a8NylOJUjESakIiYjUo9yiMsYu2MCRvBJC/VuweNK1BLTS1epFnJWKUA00R0hE7HWmpJxxCzZy4GQxwX7eLJ4UR5Cft9mxROQSNEeoFpojJCJ1UVBawbgFG9l+NJ92Pl588LsBdAloaXYsEZelOUIiIg2kpLySiW9tYvvRfNq08GDxpDiVIJFG4rLXcebk5JCTk4PNZqu2PTIy8opDiYg0FqUVVia9vZnNh0/j492Md++Po0egj9mxRKSO7C5CW7ZsISEhgT179nD+rJrFYsEwDCwWi76SQkRcRnmljQff28K6A6do6enO2xNj6dPRz+xYImIHu4vQxIkT6dGjBwsXLiQwMBCLpWl+YaC+a0xELqXSauORf//Aqn0n8fZwY+H4GK4JbWN2LBGxk92TpX18fPjhhx/o3r27ozI5FU2WFpGfs9oMEj9IZXnqcTzd3ViQ0J8be7QzO5aI/ITDJksPHTqUbdu2XVE4EZHGymYzePLjHSxPPU4zNwt/H3uNSpBII2b3qbEFCxaQkJDAzp076dOnDx4eHtXuHzVqVL2FExFxJoZh8H+f7eb9TZm4WWDOb6KJjwg0O5aIXAG7i9D69ev5/vvv+fLLLy+4T5OlRaSpMgyDF1bsZdG6QwC8dEcUt0YGmxtKRK6Y3afGpk6dyrhx4zhx4gQ2m63aj0qQiDRVc5PT+ee3BwF47vY+/LpfJ5MTiUh9sLsInTp1iscee4zAwKZ9OFhfsSEi5/3z2wO8+s1+AP54awRj4zqbnEhE6ovdRWj06NGsWrXKEVmcypQpU9i9ezebNm0yO4qImOid9YeY/eVeAJ4Y1pP7b+hiciIRqU92zxHq0aMHM2fOZO3atfTt2/eCydKPPPJIvYUTETHTB5syeXr5LgAeHtKdKUNc47IhIq7E7usIdelS87+GLBYLBw8evOJQzkTXERJxTctTj/HoklQMA+6/oQtP3dKryV5AVqQpquvnt11HhAzDYPXq1bRv357mzZtfcUgREWe0YmcWiR9swzBgbFyoSpBIE2bXHCHDMAgPD+fo0aOOyiMiYqrV+3KY+u+tWG0Go6/pyJ9/1UclSKQJs6sIubm5ER4ezqlTpxyVR0TENOsO5PK7d7dQYTW4JbIDL/46Ejc3lSCRpszuVWMvvPACTzzxBDt37nREHqeh5fMirmXL4Twmvb2Zskob8b0CmTMmmmbudv+KFJFGxu7J0m3atKGkpITKyko8PT0vmCuUl5dXrwHNpsnSIk3fjqP53PPGBgrLKhkYHsAb9/XH28Pd7FgicgUcMlkaYM6cOVeSS0TEqezNKuDeNzdSWFZJbBd/Xr9XJUjEldhdhBISEhyRQ0SkwR04WcS4BRs5U1JBdEhr3hwfQ3NPlSARV2J3EQKwWq0sW7aMPXv2ANC7d29GjRqFu7t+gYhI43DkVAlj39hIblE5vYN9eXtiLK28LutXoog0Ynb/rU9PT+fmm2/m2LFj9OzZE4DZs2cTEhLC559/Trdu3eo9pIhIfTp+5iz3LNhAVkEp4e1b8e79cfg196j9gSLS5Ni9JOKRRx6hW7duZGZmsnXrVrZu3cqRI0fo0qWLvl5DRJxeTmEpYxds5Ojps3QJaMniSXH4t/Q0O5aImMTuI0LffvstGzZswN/fv2pb27ZteeGFF7j++uvrNZyISH3KKy5n3IKNZOQW07F1cxZPiqO9r7fZsUTERHYfEfLy8qKwsPCC7UVFRXh66l9VIuKc8s9WcO/CjezPLiLQ14t/T76W4Nb6qiARV2d3Ebr11lv57W9/y8aNGzEMA8Mw2LBhAw888ACjRo1yREYRkStSVFbJ+LdS2HW8gLYtPVk86VpC27YwO5aIOAG7i9DcuXPp1q0bAwYMwNvbG29vb66//nq6d+/Oa6+95oiMptCVpUWahrPlVu5ftIkfjpzBr7kH702Ko3v7VmbHEhEnYfeVpc9LS0tj7969APTq1Yvu3bvXazBnoStLizReZZVWJr29me/ScvHxasbiyXFEdmptdiwRaQAOu7L0eeHh4YSHh1/uw0VEHKrCamPK4h/4Li2X5h7uvDUhRiVIRC5gdxGyWq0sWrSI5ORkcnJysNls1e5fuXJlvYUTEbkcVpvBo0tS+WZPNp7N3FiY0J/+Yf61P1BEXI7dRWjatGksWrSIW265hT59+mCxWByRS0TksthsBr//cDufbz+Bh7uFf47rx3XdA8yOJSJOyu4i9P777/PBBx9w8803OyKPiMhlMwyDPy7fyUdbj+LuZuFvd1/NkKvamx1LRJyY3avGPD09m+zEaBFpvAzD4NnP97B44xEsFnjlriiG9+lgdiwRcXJ2F6HHH3+c1157jctcbCYi4hCvfL2fhWszAPjL6Eh+Fd3R5EQi0hjYfWps7dq1rFq1ii+//JLevXvj4VH9iwr/85//1Fs4EZG6SFqVzt9WpgMwa1Rv7ooJMTmRiDQWdheh1q1bc/vttzsii4iI3RauzeCl/+4DYOaIq0i4LszcQCLSqNhdhN566y1H5BARsdvijYf582e7AXg0PpzfDepmciIRaWzsniMkIuIMPtpylKeW7QTgd4O6Mm2oLvAqIvZr8kUoMzOTwYMHExERQWRkJEuXLjU7kohcoc+3n+CJD7dhGJAwoDMzhl+la5qJyGW57K/YaCyaNWvGnDlziI6OJisri379+nHzzTfTsmVLs6OJyGX4Znc2097/AZsBY/qH8MzI3ipBInLZmnwR6tChAx06nLuWSFBQEAEBAeTl5akIiTRC36Wd5KHFW6m0GfwqOpjnR/fFzU0lSEQun12nxioqKhg6dChpaWn1FmDNmjWMHDmS4OBgLBYLy5Ytu2CfpKQkwsLC8Pb2Ji4ujpSUlMt6rS1btmC1WgkJ0dJakcZm48FTTH5nM+VWG8N7B/HXO6NwVwkSkStkVxHy8PBg+/bt9RqguLiYqKgokpKSLnr/kiVLSExM5JlnnmHr1q1ERUUxbNgwcnJyqvaJjo6mT58+F/wcP368ap+8vDzuu+8+Xn/99XrNLyKO98OR00xctInSChtDerZj7t1X08y9yU9xFJEGYDHsvET0Y489hpeXFy+88EL9h7FY+Pjjj7ntttuqtsXFxRETE8O8efMAsNlshISEMHXqVGbMmFGn5y0rK+MXv/gFkydP5t57761137KysqrbBQUFhISEkJ+fj6+vr/1vSkSuyK7j+dz9+gYKSiu5rltb3hwfg7eHu9mxRMTJFRQU4OfnV+vnt91zhCorK3nzzTf55ptv6Nev3wVzbV555RX709agvLycLVu2MHPmzKptbm5uxMfHs379+jo9h2EYjB8/nptuuqnWEgQwe/ZsZs2addmZRaT+pGUXcu/CFApKK+nfuQ1v3NdfJUhE6pXdRWjnzp1cc801AOzfv7/affW9ciM3Nxer1UpgYGC17YGBgezdu7dOz/H999+zZMkSIiMjq+Yfvfvuu/Tt2/ei+8+cOZPExMSq2+ePCIlIw8rILeaeBRvJKy4nspMfb06IoaVXk1/fISINzO7fKqtWrXJEDoe54YYbsNlsdd7fy8sLLy8vkpKSSEpKwmq1OjCdiFzM0dMljH1jAycLy7gqyId3Jsbi6+1R+wNFROx0RbMNjx49ytGjR+srywUCAgJwd3cnOzu72vbs7GyCgoIc9roAU6ZMYffu3WzatMmhryMi1WXll3LPGxs5nl9Kt3YteW9SHK1beJodS0SaKLuLkM1m4//+7//w8/Ojc+fOdO7cmdatW/PnP//ZriMvdeHp6Um/fv1ITk6u9vrJyckMGDCgXl9LRMyXW1TG2AUbOJJXQqh/CxZPupaAVl5mxxKRJszuU2NPPvkkCxcu5IUXXuD6668HYO3atfzpT3+itLSU5557zq7nKyoqIj09vep2RkYGqamp+Pv7ExoaSmJiIgkJCfTv35/Y2FjmzJlDcXExEyZMsDe6XXRqTKRhnSkpZ9yCjRw4WUywnzeLJ8UR5OdtdiwRaeLsXj4fHBzM/PnzGTVqVLXty5cv56GHHuLYsWN2BVi9ejVDhgy5YHtCQgKLFi0CYN68ebz00ktkZWURHR3N3LlziYuLs+t1Llddl9+JyOUrKK1g3IKNbD+aTzsfLz743QC6BOjq7yJy+er6+W13EfL29mb79u306NGj2vZ9+/YRHR3N2bNnLy+xk1IREnGskvJK7luYwubDp/Fv6cmS315LeKCP2bFEpJGr6+e33XOEoqKiqi5u+FPz5s0jKirK3qcTERdWWmFl0tub2Xz4NL7ezXhnYqxKkIg0KLvnCL344ovccsstfPPNN1UTltevX09mZiZffPFFvQc0i+YIiThWeaWNB9/bwroDp2jp6c6iibH06ehndiwRcTF2nxoDOH78OElJSVUXNezVqxcPPfQQwcHB9R7QbDo1JlL/Kq02Hv7XD6zYlYW3hxuLJsRybde2ZscSkSbEIV+xUVFRwfDhw5k/f77dq8NERACsNoPHl25jxa4sPN3deP3e/ipBImIa0799XkRch81m8OTHO1ieepxmbhb+PvYabuzRzuxYIuLC7J4sPW7cOBYuXOiILE4lKSmJiIgIYmJizI4i0iQYhsH/fbab9zdl4maBOb+JJj4isPYHiog4kN1zhKZOnco777xDeHi4w7993hlojpDIlTMMgxdW7OWf3x4E4K93RvHrfp1MTiUiTZlD5ghBw377vIg0DXOT06tK0HO391EJEhGnYVcRslqtzJo1i759+9KmTRtHZRKRJuSf3x7g1W/O/aPpj7dGMDaus8mJRET+x645Qu7u7vzyl7/kzJkzDorjPDRHSOTKvbP+ELO/PHeZjSeG9eT+G7qYnEhEpDq7J0v36dOHgwcPOiKLU5kyZQq7d+9m06ZNZkcRaZQ+2JTJ08t3AfDwkO5MGdLd5EQiIheyuwg9++yzTJ8+nc8++4wTJ05QUFBQ7UdEZHnqMf7wn3OX2rj/hi48/ssetTxCRMQcdq8ac3P7X3f66eRowzCwWCxN7isptGpMxD4rdmYx5V9bsdoMxsaF8uxtfbSQQkQanMNWja1ateqKgolI07VqXw5T/32uBI2+piN//pVKkIg4N7uL0KBBgxyRQ0QauXXpuTzw7hYqrAa3RHbgxV9H4uamEiQizs3uOUIA3333HePGjeO6667j2LFjALz77rusXbu2XsOZSavGROpu86E8Jr2zmbJKG/G9ApkzJppm7pf160VEpEHZ/Zvqo48+YtiwYTRv3pytW7dSVlYGQH5+Ps8//3y9BzSLVo2J1M32o2eY8NYmSsqtDAwPYN49V+OhEiQijcRlrRqbP38+b7zxBh4eHlXbr7/+erZu3Vqv4UTEue3NKuC+N1MoLKsktos/r9/bH28Pd7NjiYjUmd1FaN++fdx4440XbPfz83OJCy2KyDkHThYxbsFGzpRUEB3SmjfHx9DcUyVIRBoXu4tQUFAQ6enpF2xfu3YtXbt2rZdQIuLcjpwqYewbG8ktKqd3sC9vT4yllZfday9ERExndxGaPHky06ZNY+PGjVgsFo4fP87ixYuZPn06Dz74oCMyiogTOX7mLPcs2EBWQSnh7Vvx7v1x+DX3qP2BIiJOyO5/ws2YMQObzcbQoUMpKSnhxhtvxMvLi+nTpzN16lRHZBQRJ5FTWMrYBRs5evosXQJasnhSHP4tPc2OJSJy2ey+svR55eXlpKenU1RUREREBK1atarvbKZKSkoiKSkJq9XK/v37dWVpcXl5xeX85vX17M8uomPr5ix9YADBrZubHUtE5KLqemXpyy5CrkJfsSEC+WcruOeNDew6XkCgrxdLf3cdoW1bmB1LRKRGdf381sU+ROSSisoqGf9WCruOFxDQypPFk65VCRKRJkNFSERqdLbcyv2LNvHDkTO0buHBu/fH0b190zoNLiKuTUVIRC6qrNLKb9/dzMaMPHy8mvHOxFh6ddDpYRFpWlSEROQCFVYbUxb/wHdpuTT3cOetCTFEdmptdiwRkXqnIiQi1VhtBo8uSeWbPdl4NnNjYUJ/+of5mx1LRMQhVIREpIrNZvD7D7fz+fYTeLhb+Oe4flzXPcDsWCIiDqMiJCIAGIbBH5fv5KOtR3F3s/C3u69myFXtzY4lIuJQKkIigmEYPPv5HhZvPILFAq/cFcXwPh3MjiUi4nAqQjVISkoiIiKCmJgYs6OIONwrX+9n4doMAP4yOpJfRXc0OZGISMPQlaVroStLS1OXtCqdl/67D4BZo3qTcF2YuYFEROqBriwtIrVauDajqgTNHHGVSpCIuBwVIREXtXjjYf782W4AHo0P53eDupmcSESk4akIibigj7Yc5allOwH43aCuTBsabnIiERFzqAiJuJjPt5/giQ+3YRiQMKAzM4ZfhcViMTuWiIgpVIREXMg3u7OZ9v4P2AwY0z+EZ0b2VgkSEZemIiTiIr5LO8lDi7dSaTP4VXQwz4/ui5ubSpCIuDYVIREXsPHgKSa/s5lyq43hvYP4651RuKsEiYioCIk0dT8cOc3ERZsorbAxpGc75t59Nc3c9VdfRASgmdkBRKT+WG0GKRl55BSW0t7Hmxae7iS8mUJxuZXrurXlH+P64dlMJUhE5DwVIZEmYsXOE8z6dDcn8kurtlksYBjQv3Mb3rivP94e7iYmFBFxPk3+n4Znzpyhf//+REdH06dPH9544w2zI4nUuxU7T/Dge1urlSA4V4IA7okLpaWX/t0jIvJzTf43o4+PD2vWrKFFixYUFxfTp08fRo8eTdu2bc2OJlIvrDaDWZ/upqYvDbQAL/13H7+K7qgJ0iIiP9Pkjwi5u7vTokULAMrKyjAMA33PrDQlKRl5FxwJ+ikDOJFfSkpGXsOFEhFpJEwvQmvWrGHkyJEEBwdjsVhYtmzZBfskJSURFhaGt7c3cXFxpKSk2PUaZ86cISoqik6dOvHEE08QEBBQT+lFzJdTWHMJupz9RERcielFqLi4mKioKJKSki56/5IlS0hMTOSZZ55h69atREVFMWzYMHJycqr2OT//5+c/x48fB6B169Zs27aNjIwM/vWvf5Gdnd0g702kIXjVcRVYex9vBycREWl8LIYTnSeyWCx8/PHH3HbbbVXb4uLiiImJYd68eQDYbDZCQkKYOnUqM2bMsPs1HnroIW666SbuuOOOi95fVlZGWVlZ1e2CggJCQkLIz8/H19fX7tcTcaSvd2fzhw+3kVdSUeM+FiDIz5u1f7hJc4RExGUUFBTg5+dX6+e36UeELqW8vJwtW7YQHx9ftc3NzY34+HjWr19fp+fIzs6msLAQgPz8fNasWUPPnj1r3H/27Nn4+flV/YSEhFzZmxBxgOKySmZ8tJ3J72wmr6SCjq3PHe35ec05f/uZkREqQSIiF+HURSg3Nxer1UpgYGC17YGBgWRlZdXpOQ4fPszAgQOJiopi4MCBTJ06lb59+9a4/8yZM8nPz6/6yczMvKL3IFLfthw+zc1zv+P9TZlYLPDbG7uycvpg5o+7hiC/6qe/gvy8+ce4axjep4NJaUVEnFuTXz4fGxtLampqnff38vLCy8vLcYFELlOF1cbc5DSSVqVjMyDYz5u/3hXNgG7nLgUxvE8HfhERVO3K0rFd/HUkSETkEpy6CAUEBODu7n7B5Obs7GyCgoIc+tpJSUkkJSVhtVod+joidXHgZBGPLUll+9F8AG6/uiN/GtUbv+Ye1fZzd7NUFSMREamdU58a8/T0pF+/fiQnJ1dts9lsJCcnM2DAAIe+9pQpU9i9ezebNm1y6OuIXIphGLyz/hC3zP2O7Ufz8Wvuwbx7rubVMdEXlCAREbGf6UeEioqKSE9Pr7qdkZFBamoq/v7+hIaGkpiYSEJCAv379yc2NpY5c+ZQXFzMhAkTTEwt4ng5BaU88eF2vt1/EoCB4QG8dEfUBfOARETk8plehDZv3syQIUOqbicmJgKQkJDAokWLGDNmDCdPnuTpp58mKyuL6OhoVqxYccEE6vqmU2NiphU7TzDzPzs4XVKBVzM3Zoy4ioQBYbhpvo+ISL1yqusIOaO6XodApD4Ullbwp09289HWowD0DvZlzphowgN9TE4mItK41PXz2/QjQiJyTkpGHokfpHL09FncLPDAoG48Gt8DzzpeOVpEROynIlQDnRqThlJeaePVb/Yz/9sDGAZ0atOcV8dEExPmb3Y0EZEmT6fGaqFTY+JI+7MLefT9VHafKADgzn6deHpkBD7eWhEmInIldGpMxInZbAaL1h3ihRV7Ka+00aaFB7NH99UVoEVEGpiKkEgDO5F/lieWbmdtei4Ag3u248VfR9LeV8viRUQamopQDTRHSBzh023HefLjHRSUVuLt4caTt0QwLi4Ui0XL4kVEzKA5QrXQHCGpD/lnK3h6+U6Wpx4HIKqTH6+MiaZbu1YmJxMRaZo0R0jESaw7kMv0D7ZxPL8UdzcLU4Z0Z+pN3fFw17J4ERGzqQiJOEhphZW/frWPBWszMAzo3LYFr46J5prQNmZHExGRH6kIiTjAnhMFPLYklb1ZhQDcHRvKU7f0oqWX/sqJiDgT/VaugSZLy+Ww2gwWfHeQv361n3KrjYBWnrwwOpL4CMd+N56IiFweTZauhSZLS10dPV3C4x9sY2NGHgDxvdrzwq8jCWjlZXIyERHXo8nSIg3EMAyWpR7j6WW7KCyrpIWnO0/fGsGYmBAtixcRcXIqQiJX4ExJOU9+vJPPd5wA4OrQ1rx6VzRhAS1NTiYiInWhIiRymb5LO8n0pdvILijD3c3Co0PDeXBwN5ppWbyISKOhIlQDTZaWmpRWWHnhy70sWncIgK7tWvLqXdFEhbQ2NZeIiNhPk6VrocnS8lM7j+Xz6JJU0nOKALhvQGdmjuhFc093k5OJiMhPabK0SD2y2gzmf3uAV7/eT6XNoJ2PFy/dEcngnu3NjiYiIldARUikFpl5JTy2JJXNh08DMLx3EM+P7ot/S0+Tk4mIyJVSERKpgWEYLN1ylFmf7KK43Eorr2b8aVRvfn1NRy2LFxFpIlSERC7iVFEZ/+/jHfx3VzYAMWFteOWuaEL8W5icTERE6pOKkMjPrNqbwxMfbie3qAwPdwuP/aIHv7uxG+5uOgokItLUqAiJ/KikvJLnv9jDexuOABDevhWvjommT0c/k5OJiIijqAjVQNcRci2pmWdIXJLKwdxiACZe34XfD++Jt4eWxYuINGW6jlAtdB2hpq3SaiNp1QHmrkzDajMI8vXm5TujuCE8wOxoIiJyBXQdIZFaZOQW89iSVFIzzwBwa2QHnr2tD61baFm8iIirUBESl2MYBv9OyeTPn+3mbIUVH+9mPHtbH0ZFBWtZvIiIi1EREpdysrCMGR9tJ3lvDgADurbl5bui6Ni6ucnJRETEDCpC4jK+3p3NjI+2c6q4HE93N54Y1pP7b+iCm5bFi4i4LBUhafKKyyr5v093s2RzJgBXBfkw5zfRXBWkye8iIq5ORUiatC2HT/PYklSO5JVgscDkgV15/Jc98GqmZfEiIqIiJE1UhdXG3OQ0klalYzMg2M+bv94VzYBubc2OJiIiTkRFSJqc9JwiHluSyo5j+QDcfnVH/jSqN37NPUxOJiIizkZFqAa6snTjYxgG7244zPNf7KG0woZfcw+eu70Pt0YGmx1NRESclK4sXQtdWbpxyC4o5YkPt7Nm/0kABoYH8NIdUQT5eZucTEREzKArS4vL+HLHCWZ+vIMzJRV4NXNjxoirSBgQpmXxIiJSKxUhabQKSyv40ye7+WjrUQB6B/syZ0w04YE+JicTEZHGQkVIGqWUjDweW5LKsTNncbPAA4O68Wh8DzybuZkdTUREGhEVIWlUyiqtvPp1Gv9ccwDDgE5tmvPqmGhiwvzNjiYiIo2QipA0GvuzC5n2fip7ThQAcGe/Tjw9MgIfby2LFxGRy6MiJE7PZjN4a90h/rJiL+WVNtq08GD26L4M79PB7GgiItLIqQiJUzuRf5bpS7fxffopAAb3bMeLv46kva+WxYuIyJVTERKn9cm24zz18Q4KSivx9nDjyVsiGBcXisWiZfEiIlI/VITE6eSfreDp5TtZnnocgKhOfrwyJppu7VqZnExERJoaFSFxKuvSc3l86TZO5Jfi7mZhypDuTL2pOx7uWhYvIiL1z2U+XUpKSujcuTPTp083O4pcRGmFlWc/2809CzZyIr+Uzm1bsPSBAST+oodKkIiIOIzLHBF67rnnuPbaa82OIRex+3gBjy1JZV92IQB3x4by1C29aOnlMn88RUTEJC7xSZOWlsbevXsZOXIkO3fuNDuO/MhqM1jw3UH++tV+yq02Alp58sLoSOIjAs2OJiIiLsL0cw5r1qxh5MiRBAcHY7FYWLZs2QX7JCUlERYWhre3N3FxcaSkpNj1GtOnT2f27Nn1lFjqw9HTJdzzxgZmf7mXcquN+F7tWfHojSpBIiLSoEw/IlRcXExUVBQTJ05k9OjRF9y/ZMkSEhMTmT9/PnFxccyZM4dhw4axb98+2rdvD0B0dDSVlZUXPParr75i06ZN9OjRgx49erBu3bpa85SVlVFWVlZ1u6Cg4ArenfycYRh8/MMxnlm+i8KySlp4uvP0rRGMiQnRsngREWlwFsMwDLNDnGexWPj444+57bbbqrbFxcURExPDvHnzALDZbISEhDB16lRmzJhR63POnDmT9957D3d3d4qKiqioqODxxx/n6aefvuj+f/rTn5g1a9YF2/Pz8/H19b28NyYAnCkp58mPd/L5jhMAXB3amlfviiYsoKXJyUREpKkpKCjAz8+v1s9vpy5C5eXltGjRgg8//LBaOUpISODMmTMsX77crudftGgRO3fu5OWXX65xn4sdEQoJCVERukLfpZ1k+tJtZBeU0czNwrSh4Tw4uBvNtCJMREQcoK5FyPRTY5eSm5uL1WolMLD6vJHAwED27t3rkNf08vLCy8vLIc/tikorrLzw5V4WrTsEQNd2LZkzJprITq1NzSUiIgJOXoTq2/jx4+u8b1JSEklJSVitVscFauJ2Hstn2vs/cOBkMQD3DejMzBG9aO7pbnIyERGRc5y6CAUEBODu7k52dna17dnZ2QQFBTn0tadMmcKUKVOqDq1J3VltBvO/PcCrX++n0mbQzseLl+6IZHDP9mZHExERqcapJ2h4enrSr18/kpOTq7bZbDaSk5MZMGCAicmkJkdOlTDmn+t56b/7qLQZDO8dxH8fvVElSEREnJLpR4SKiopIT0+vup2RkUFqair+/v6EhoaSmJhIQkIC/fv3JzY2ljlz5lBcXMyECRMcmkunxuxjGAZLNx9l1qe7KC630sqrGX8a1ZtfX9NRy+JFRMRpmb5qbPXq1QwZMuSC7QkJCSxatAiAefPm8dJLL5GVlUV0dDRz584lLi6uQfLVdda5KztVVMbM/+zgq93nTmHGhLXhlbuiCfFvYXIyERFxVY1y+bwzUhG6tFV7c3jiw+3kFpXh4W7hsV/04Hc3dsPdTUeBRETEPE1i+byZdGrs0krKK3nu8z0s3ngEgPD2rXh1TDR9OmpiuYiINB46IlQLHRG6UGrmGR5bkkpG7rll8ROv78Lvh/fE20PL4kVExDnoiJDUu0qrjXmr0vnbynSsNoMgX29evjOKG8IDzI4mIiJyWVSEpE4ycot5bEkqqZlnALg1sgPP3taH1i08zQ0mIiJyBVSEaqA5QucYhsG/Uo7w7Gd7OFthxce7Gc/e1odRUcFaFi8iIo2e5gjVwpXnCJ0sLGPGR9tJ3psDwICubXn5rig6tm5ucjIREZFL0xwhuSJf7cpixn92kFdcjqe7G08M68n9N3TBTcviRUSkCVERkmqKyir586e7WbI5E4CrgnyY85torgpyraNhIiLiGlSEpMqWw3k8tmQbR/JKsFhg8sCuPP7LHng107J4ERFpmlSEauBKk6UrrDZe+yaNv69Ox2ZAx9bNefnOKAZ0a2t2NBEREYfSZOlaNPXJ0uk5RTy2JJUdx/IBuP3qjsz6VW98vT1MTiYiInL5NFlaLskwDN5Zf5jnv9hDWaUNv+YePHd7H26NDDY7moiISINREXJB2QWlPPHhdtbsPwnAwPAAXrojiiA/b5OTiYiINCwVIRfzxY4T/L+Pd3CmpAKvZm7MGHEVCQPCtCxeRERckoqQiygoreBPn+ziP1uPAdA72Jc5Y6IJD/QxOZmIiIh5VIRq0JRWjW08eIrED7Zx7MxZ3CzwwKBuPBrfA89mbmZHExERMZVWjdWiMa8aK6u08srX+3l9zUEMAzq1ac6rY6KJCfM3O5qIiIhDadWYi9ufXci091PZc6IAgDv7deLpkRH4aFm8iIhIFRWhJsZmM3hr3SH+smIv5ZU22rTwYPbovgzv08HsaCIiIk5HRagJOZF/lulLt/F9+ikABvdsx4t3RNLeR8viRURELkZFqIn4ZNtxnvp4BwWllXh7uPHkLRGMiwvFYtGyeBERkZqoCDVy+SUVPP3JTpanHgcgqpMfr4yJplu7ViYnExERcX4qQjVoDMvn16Xn8vjSbZzIL8XdzcKUId2ZelN3PNy1LF5ERKQutHy+Fs64fL60wspL/93HwrUZAHRu24JXx0RzTWgbk5OJiIg4By2fb6J2Hy/gsSWp7MsuBODu2FCeuqUXLb30n1JERMRe+vRsJKw2gwXfHeSvX+2n3GojoJUnL4yOJD4i0OxoIiIijZaKUCNw9HQJiR9sIyUjD4D4XoG88Ou+BLTyMjmZiIhI46Yi5MQMw+A/W4/xp092UVhWSQtPd56+NYIxMSFaFi8iIlIPVISc1Onicp5ctoMvdmQBcHVoa169K5qwgJYmJxMREWk6VISc0Jr9J5m+dBs5hWU0c7MwbWg4Dw7uRjMtixcREalXKkImsNoMUjLyyCkspb2PN7Fd/HF3s3C23MoLX+7h7fWHAejariVzxkQT2am1uYFFRESaKBWhBrZi5wlmfbqbE/mlVds6+Hkz4bowlmzO5MDJYgDuG9CZmSN60dzT3ayoIiIiTZ6KUA0ccWXpFTtP8OB7W/n5FSxP5Jfy/Jd7AWjn48VLd0QyuGf7entdERERuThdWboW9XVlaavN4Ia/rKx2JOjnvJu58d0fbqKdj5bFi4iIXIm6fn5r9m0DScnIu2QJAiittJGeU9RAiURERERFqIHkFF66BNm7n4iIiFw5FaEG0t7Hu173ExERkSunItRAYrv408HPm5quB23h3Oqx2C7+DRlLRETEpakINRB3NwvPjIwAuKAMnb/9zMgI3N301RkiIiINRUWoAQ3v04F/jLuGIL/qp7+C/Lz5x7hrGN6ng0nJREREXJOuI9TAhvfpwC8igi56ZWkRERFpWCpCJnB3szCgW1uzY4iIiLg8nRoTERERl6UiJCIiIi7LJU6NhYWF4evri5ubG23atGHVqlVmRxIREREn4BJFCGDdunW0atXK7BgiIiLiRHRqTERERFyW6UVozZo1jBw5kuDgYCwWC8uWLbtgn6SkJMLCwvD29iYuLo6UlBS7XsNisTBo0CBiYmJYvHhxPSUXERGRxs70U2PFxcVERUUxceJERo8efcH9S5YsITExkfnz5xMXF8ecOXMYNmwY+/bto3379gBER0dTWVl5wWO/+uorgoODWbt2LR07duTEiRPEx8fTt29fIiMjHf7eRERExLlZDMMwzA5xnsVi4eOPP+a2226r2hYXF0dMTAzz5s0DwGazERISwtSpU5kxY4bdr/HEE0/Qu3dvxo8ff9H7y8rKKCsrq7pdUFBASEgI+fn5+Pr62v16IiIi0vAKCgrw8/Or9fPb9FNjl1JeXs6WLVuIj4+v2ubm5kZ8fDzr16+v03MUFxdTWFgIQFFREStXrqR379417j979mz8/PyqfkJCQq7sTYiIiIjTMv3U2KXk5uZitVoJDAystj0wMJC9e/fW6Tmys7O5/fbbAbBarUyePJmYmJga9585cyaJiYlVt/Pz8wkNDaWgoOAy3oGIiIiY4fzndm0nvpy6CNWHrl27sm3btjrv7+XlhZeXV9Xt8wOpI0MiIiKNT2FhIX5+fjXe79RFKCAgAHd3d7Kzs6ttz87OJigoqEEyBAcHk5mZiY+PDxZL9S9GjYmJYdOmTZfcVtPt83OPMjMzHTL36GLZ6usxl9qvpvs0VnW/70rGCnDoeGms6u5yxqquj3PUWP18W0ON1aVyX+ljGmqsfnq7KY7Vpe531t/vhmFQWFhIcHDwJfdz6iLk6elJv379SE5OrppAbbPZSE5O5uGHH26QDG5ubnTq1Omi97m7u1/wH+7n22q77evr65C/KBfLVl+PudR+Nd2nsar7ffUxVuCY8dJY1d3ljFVdH+eosfr5toYaq5peqz4e01BjdbHbTWmsLnW/M/9+v9SRoPNML0JFRUWkp6dX3c7IyCA1NRV/f39CQ0NJTEwkISGB/v37Exsby5w5cyguLmbChAkmpj5nypQptW6r7bajXM7r1PUxl9qvpvs0VnW/T2NV9/ua2ljV9XGOGqufb2uosbrc13Kmsaprnvpgxlhd6n5n/ntYF6Yvn1+9ejVDhgy5YHtCQgKLFi0CYN68ebz00ktkZWURHR3N3LlziYuLa+Ck9auuy/pEY2UvjVfdaazqTmNVdxqrunOGsTL9iNDgwYNrndH98MMPN9ipsIbi5eXFM888U21itlycxso+Gq+601jVncaq7jRWdecMY2X6ESERERERszj1BRVFREREHElFSERERFyWipCIiIi4LBUhERERcVkqQiIiIuKyVIQaiVdffZXevXsTERHBI488UuslB1zVvn37iI6Orvpp3rw5y5YtMzuW08rIyGDIkCFERETQt29fiouLzY7ktMLCwoiMjCQ6Ovqi1z6T6kpKSujcuTPTp083O4rTOnPmDP379yc6Opo+ffrwxhtvmB3JqWVmZjJ48GAiIiKIjIxk6dKl9fK8Wj7fCJw8eZJrr72WXbt24eHhwY033sjLL7/MgAEDzI7m1IqKiggLC+Pw4cO0bNnS7DhOadCgQTz77LMMHDiQvLw8fH19adbM9MuLOaWwsDB27txJq1atzI7SKDz55JOkp6cTEhLCyy+/bHYcp2S1WikrK6NFixYUFxfTp08fNm/eTNu2bc2O5pROnDhBdnY20dHRZGVl0a9fP/bv33/Fv991RKiRqKyspLS0lIqKCioqKmjfvr3ZkZzeJ598wtChQ1WCanC+WA8cOBAAf39/lSCpF2lpaezdu5cRI0aYHcWpubu706JFCwDKysowDENH+y+hQ4cOREdHAxAUFERAQAB5eXlX/LwqQvVgzZo1jBw5kuDgYCwWy0VPxSQlJREWFoa3tzdxcXGkpKTU+fnbtWvH9OnTCQ0NJTg4mPj4eLp161aP76DhOHqsfuqDDz5gzJgxV5jYPI4eq7S0NFq1asXIkSO55ppreP755+sxfcNqiD9XFouFQYMGERMTw+LFi+specNriLGaPn06s2fPrqfE5mmIsTpz5gxRUVF06tSJJ554goCAgHpK3/Aa8vf7li1bsFqthISEXGFqFaF6UVxcTFRUFElJSRe9f8mSJSQmJvLMM8+wdetWoqKiGDZsGDk5OVX7nD9H/POf48ePc/r0aT777DMOHTrEsWPHWLduHWvWrGmot1evHD1W5xUUFLBu3Tpuvvlmh78nR3H0WFVWVvLdd9/x97//nfXr1/P111/z9ddfN9Tbq1cN8edq7dq1bNmyhU8++YTnn3+e7du3N8h7q2+OHqvly5fTo0cPevTo0VBvyWEa4s9V69at2bZtGxkZGfzrX/8iOzu7Qd6bIzTU7/e8vDzuu+8+Xn/99foJbki9AoyPP/642rbY2FhjypQpVbetVqsRHBxszJ49u07P+cEHHxgPPfRQ1e0XX3zR+Mtf/lIvec3kiLE675133jHGjh1bHzGdgiPGat26dcYvf/nLqtsvvvii8eKLL9ZLXjM58s/VedOnTzfeeuutK0jpHBwxVjNmzDA6depkdO7c2Wjbtq3h6+trzJo1qz5jm6Ih/lw9+OCDxtKlS68kptNw1HiVlpYaAwcONN555536imroiJCDlZeXs2XLFuLj46u2ubm5ER8fz/r16+v0HCEhIaxbt47S0lKsViurV6+mZ8+ejopsmvoYq/Ma+2mx2tTHWMXExJCTk8Pp06ex2WysWbOGXr16OSqyaepjrIqLiyksLATOTcJfuXIlvXv3dkheM9XHWM2ePZvMzEwOHTrEyy+/zOTJk3n66acdFdk09TFW2dnZVX+u8vPzWbNmTZP83Q71M16GYTB+/Hhuuukm7r333nrLpiLkYLm5uVitVgIDA6ttDwwMJCsrq07Pce2113LzzTdz9dVXExkZSbdu3Rg1apQj4pqqPsYKzv1CSUlJYdiwYfUd0WnUx1g1a9aM559/nhtvvJHIyEjCw8O59dZbHRHXVPUxVtnZ2dxwww1ERUVx7bXXct999xETE+OIuKaqr7+DrqA+xurw4cMMHDiQqKgoBg4cyNSpU+nbt68j4pquPsbr+++/Z8mSJSxbtqzqEik7duy44mxaItJIPPfcczz33HNmx2gU/Pz8GvV59oY0YsQIreypg65du7Jt2zazYzQ648ePNzuCU4uNjSU1NdXsGI3GDTfcgM1mq/fn1REhBwsICMDd3f2CD+bs7GyCgoJMSuWcNFZ1p7GqO41V3Wms6k5jZR9nHi8VIQfz9PSkX79+JCcnV22z2WwkJyfrgog/o7GqO41V3Wms6k5jVXcaK/s483jp1Fg9KCoqIj09vep2RkYGqamp+Pv7ExoaSmJiIgkJCfTv35/Y2FjmzJlDcXExEyZMMDG1OTRWdaexqjuNVd1prOpOY2WfRjte9bb+zIWtWrXKAC74SUhIqNrnb3/7mxEaGmp4enoasbGxxoYNG8wLbCKNVd1prOpOY1V3Gqu601jZp7GOl75rTERERFyW5giJiIiIy1IREhEREZelIiQiIiIuS0VIREREXJaKkIiIiLgsFSERERFxWSpCIiIi4rJUhERERMRlqQiJSKMzePBgHn30UbNjiEgToCIkIqZSqRERM6kIiYjLsVqt2Gw2s2OIiBNQERIR04wfP55vv/2W1157DYvFgsVi4dChQ3z77bfExsbi5eVFhw4dmDFjBpWVlTU+T1lZGdOnT6djx460bNmSuLg4Vq9eXXX/okWLaN26NZ988gkRERF4eXlx5MgRNm3axC9+8QsCAgLw8/Nj0KBBbN26tdpzWywWFixYwO23306LFi0IDw/nk08+qbbPrl27uPXWW/H19cXHx4eBAwdy4MCBqvsXLFhAr1698Pb25qqrruLvf/97/QygiFwxFSERMc1rr73GgAEDmDx5MidOnODEiRN4eHhw8803ExMTw7Zt2/jHP/7BwoULefbZZ2t8nocffpj169fz/vvvs337du68806GDx9OWlpa1T4lJSX85S9/YcGCBezatYv27dtTWFhIQkICa9euZcOGDYSHh3PzzTdTWFhY7flnzZrFXXfdxfbt27n55psZO3YseXl5ABw7dowbb7wRLy8vVq5cyZYtW5g4cWJVcVu8eDFPP/00zz33HHv27OH555/nj3/8I2+//bYDRlRE7GbeF9+LiBjGoEGDjGnTplXd/n//7/8ZPXv2NGw2W9W2pKQko1WrVobVar3gMYcPHzbc3d2NY8eOVXveoUOHGjNnzjQMwzDeeustAzBSU1MvmcVqtRo+Pj7Gp59+WrUNMJ566qmq20VFRQZgfPnll4ZhGMbMmTONLl26GOXl5Rd9zm7duhn/+te/qm3785//bAwYMOCSWUSkYTQzuYeJiFSzZ88eBgwYgMViqdp2/fXXU1RUxNGjRwkNDa22/44dO7BarfTo0aPa9rKyMtq2bVt129PTk8jIyGr7ZGdn89RTT7F69WpycnKwWq2UlJRw5MiRavv99HEtW7bE19eXnJwcAFJTUxk4cCAeHh4XvJfi4mIOHDjA/fffz+TJk6u2V1ZW4ufnV9chEREHUhESkUatqKgId3d3tmzZgru7e7X7WrVqVfX/mzdvXq1cASQkJHDq1Clee+01OnfujJeXFwMGDKC8vLzafj8vORaLpWqydfPmzS+ZDeCNN94gLi6u2n0/zyoi5lAREhFTeXp6YrVaq2736tWLjz76CMMwqorL999/j4+PD506dbrg8VdffTVWq5WcnBwGDhxo12t///33/P3vf+fmm28GIDMzk9zcXLueIzIykrfffpuKiooLClNgYCDBwcEcPHiQsWPH2vW8ItIwNFlaREwVFhbGxo0bOXToELm5uTz00ENkZmYydepU9u7dy/Lly3nmmWdITEzEze3CX1k9evRg7Nix3HffffznP/8hIyODlJQUZs+ezeeff37J1w4PD+fdd99lz549bNy4kbFjx17yCM/FPPzwwxQUFPCb3/yGzZs3k5aWxrvvvsu+ffuAcxOtZ8+ezdy5c9m/fz87duzgrbfe4pVXXrHrdUTEMVSERMRU06dPx93dnYiICNq1a0dFRQVffPEFKSkpREVF8cADD3D//ffz1FNP1fgcb731Fvfddx+PP/44PXv25LbbbmPTpk0XzCf6uYULF3L69GmuueYa7r33Xh555BHat29vV/62bduycuVKioqKGDRoEP369eONN96oOjo0adIkFixYwFtvvUXfvn0ZNGgQixYtokuXLna9jog4hsUwDMPsECIiIiJm0BEhERERcVkqQiIiIuKyVIRERETEZakIiYiIiMtSERIRERGXpSIkIiIiLktFSERERFyWipCIiIi4LBUhERERcVkqQiIiIuKyVIRERETEZakIiYiIiMv6/6ULblXr1wz7AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "models = pybamm.lithium_ion.DFN()\n",
    "high_tol = 1e-10\n",
    "goal_rtol = 1e-4\n",
    "goal_atol = 1e-6\n",
    "t_eval = np.linspace(0, 3600, 100)\n",
    "solver = pybamm.IDAKLUSolver(atol=high_tol, rtol=high_tol)\n",
    "sim = pybamm.Simulation(models, solver=solver)\n",
    "high_tol_sol = sim.solve([0, 3600])[\"Voltage [V]\"](t_eval)\n",
    "\n",
    "tols = [1e-2, 1e-4, 1e-6, 1e-8]\n",
    "results = np.zeros(len(tols))\n",
    "for itol, tol in enumerate(tols):\n",
    "    solver = pybamm.IDAKLUSolver(atol=tol, rtol=tol)\n",
    "    sim = pybamm.Simulation(model, solver=solver)\n",
    "    sol = sim.solve([0, 3600])[\"Voltage [V]\"](t_eval)\n",
    "    results[itol] = np.sqrt(\n",
    "        np.sum(\n",
    "            ((sol - high_tol_sol) / (goal_rtol * np.abs(high_tol_sol) + goal_atol)) ** 2\n",
    "        )\n",
    "        / len(sol)\n",
    "    )\n",
    "\n",
    "# plot results on a log scale\n",
    "plt.figure()\n",
    "plt.plot(tols, results, \"o-\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xscale(\"log\")\n",
    "plt.xlabel(\"tolerance\")\n",
    "plt.ylabel(\"error norm\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us some justification that a tolerance of $10^{-4}$ is a good compromise between accuracy and computational time for this problem. However, the choice of tolerances will be highly problem-dependent, and your measure of error may be different to the one used here, so it is always best to determine the appropriate error measure and convergence study for your particular problem."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
